Final Report: 


Acoustic Fatigue Life Prediction for Nonlinear Structures 


with Multiple Resonant Modes 
NASA Grant #NAG-l-978 
April 1989 - March 1992 


//} a ^ ^ £’7 

/AJ 


R. N. Miles 

Department of Mechanical and Industrial Engineering 
State University of New York 


P. 0. 6000 

Binghamton, NY 13902-6000 


J 6 V S' ^ O 


Abstract 

This report documents an effort to develop practical and accurate methods of es- 
timating the fatigue lives of complex aerospace structures subjected to intense random 
excitations. The emphasis of the current program is to construct analytical schemes of 
performing fatigue life estimates for structures that exhibit nonlinear vibration behavior 
and that have numerous resonant modes contributing to the response. 
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I. Introduction 


This report describes results obtained over a three year period under funding from the 
Structural Acoustics, Branch at NASA Langley Research Center. The goal of this effort has 
been to develop methods for predicting the acoustic fatigue lives of aerospace structures. 
An attempt has been made to focus on certain technology areas that this investigator 
believes are in most need of further development. The two main issues that have been 
addressed are 1) how to properly account for the effects of broadband random loading in 
acoustic fatigue studies and 2) how to include nonlinear effects due to intense excitations 
in acoustic fatigue predictions. 

In acoustic fatigue studies it is important to be able to account for very complicated 
stress and strain time histories because acoustic loads often contain components over a 
broad range of frequencies. As a result, a large number of resonant modes can be excited. 
As will be discussed in the following, although acoustic fatigue problems are often charac- 
terized by extremely broadband loads, (in contrast to most fatigue problems) it is common 
practice to idealize structures as having a single degree of freedom. In this case, the re- 
sponse is narrowband. This greatly simplifies fatigue predictions but it is not consistent 
with observations. 


To aid in the development of ‘consistent 5 fatigue prediction procedures, several sections 
are devoted to the problem of identifying the proper damage cycle counting method for 
acoustic fatigue predictions. The Rainflow cycle counting method is discussed in detail 
and methods are developed to aid in its implementation. 


Nonlinear effects are often important in acoustic fatigue studies because it is the high 
level excitations that cause the most damage. Nonlinear structural response can have a 
significant impact on acoustic fatigue life in service as well as in tests. Acoustic fatigue 
testing is often performed using artificially high excitation levels in order to accelerate the 
rate of damage accumulation. Accelerated testing is necessary when trying to determine 
the fatigue characteristics of a structure that is expected to be useful for 20 years. When 
nonlinear behavior contributes to the structural response, the fatigue mechanisms can be 
modified significantly. 


In order to account for nonlinearities in fatigue predictions it is necessary to obtain 
a time domain simulation of the response. In a nonlinear structure with a large number 
of resonant modes, as often considered in acoustic fatigue studies, it is not computation- 
ally practical to perform a detailed time domain simulation of the random response. In 
this report, a numerical procedure is developed which permits a fatigue prediction for a 
nonlinear structure to be performed with roughly the same numerical effort as for a linear 
structure having the same number of degrees of freedom. The procedure is applied to 
nonlinear beams and plates and fatigue predictions are found to agree closely with results 
of more computationally intense conventional methods. 


\ 
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II. 1 Basics of Fatigue Life Estimation - Peak Counting 


For the present study we will assume that the relation between the stress am plitude, S, 
and the number of cycles to failure N, may be approximated by the well known relation, 



(II.l.l) 


where c and b are experimentally obtained constants for a given material. For the present 
discussion it will be assumed that the number of fatigue cycles experienced by the structure 
is equal to the number of positive stress peaks, or stress reversals, that occur over time. 
This approach has been used in a number of acoustic fatigue studies. An alternative cycle 
counting method, Rainflow cycle counting, will be described in later sections and will be 
compared to the older peak counting scheme. 


The assumption that damage occurs at each positive stress peak may be conveniently 
combined with the Palmgren- Miner linear damage accumulation rule [1] to produce fatigue 
life estimates. In this theory, the damage, X?,, caused by stress reversals at the stress level 
Si is 


D,(SiT= 


n(S ,) 

N(SiY 


(II.1.2) 


where n(5,) is the number of stress reversals experienced by the structure at the stress 
level Si and N(Si ) is the number of reversals required to cause failure at this stress level. 
The total damage, D m , will be the sum of the damage at all stress levels that occur, 


D m = £A(Si) = £ 


n(S j) 

mSiY 


(II.1.3) 


From equation (II.l.l), 

N (Si) = (II.1.4) 

so that equation (II. 1.3) becomes 

Dm = ~ ^r^n(Si)£j. (II. 1.5) 

t 


The form of equation (II.l.l) assumes that all stress peaks occur at positive stress 
levels as would be the case for sinusoidal loading. When a resonant system is subjected 
to random loading, however, we must allow for the possibility of stress peaks at negative 
stresses which make positive contributions to the accumulated damage. To account for this, 
we will take the absolute value of the stress level, Si , so that equation (II. 1.5) becomes 

D m = -y>(S0|S,-| 6 . (II.1.6) 

C 


Failure is predicted to occur when D m = 1. 
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Since the response and rate of damage accumulation in the structure are assumed to 
be random, the fatigue life may be estimated from the expected value of the rate of damage 
accumulation. If m(S) is the rate of occurrence of a peak with level S', and if p pea k(S ) is 
the probability density for response peaks, then the expected value of the damage rate, 
A(S), may be written as 


/ oo I C| 6 

m(5)i-Lp peoJfc (5)dS'. 

-oo C 


(II. 1.7) 


The assumption that the response process is stationary gives the mean fatigue life as 

1 


T = 


E[A]' 


(II.1.8) 


Prom equations (II. 1.7) and (II. 1.8) the problem of estimating fatigue life in a ran- 
dom system is one of properly determining the rate at which peaks occur, m(S), and the 
peak probability density, p pea k(S). The goal of the present study is to develop practical 
and accurate methods of applying equations (II. 1.7) and (II.1.8) to structures with multi- 
ple resonant modes and where the response levels are sufficiently high to elicit nonlinear 
response. 

To illustrate the application of equations (II.l. 7) and (II.1.8) consider the problem 
of predicting the fatigue life of a linear system having one resonant frequency and that is 
driven with Gaussian white noise, 


x + u 2 0 x + 0x = f(t). (II.1.9) 

In a linear system the stress is linearly related to the displacement, x by S = Kx, where 
the constant, K depends on material and geometrical properties. It may be shown that 
for Gaussian input, the response of the linear system will be Gaussian. It is also known 
that for a Gaussian random process, the peak probability density depends on the response 
spectrum. 

Assume that the response of a one degree of freedom oscillator is a narrowband process. 
Then the probability density for peaks in the time domain response is the Rayleigh density, 

Ppeak(S) = (II.1.10) 

a s 

where a 2 s is the mean square stress. Also assume that the rate at which peaks occur is a 
constant and equal to the natural frequency of the oscillator, 

( 11 - 1 * 11 ) 

Substituting equations (II.l. 10) and (II. 1.11) into equation (II. 1.7) and carrying out 
the integration gives 

E[A] = ^(v/2<r s )‘ r(t±l), (II.l. 12) 
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where the Gamma function is defined by 



for y > 0. 


Equation (II.1.12) was first obtained by Miles [2]. 


(II.l. 13) 


II. 2. Linear Systems with Multiple Degrees of Freedom 

In the case of a single degree of freedom system as discussed in the previous section, it 
is relatively straight-forward to estimate the rate at which damaging events occur in the 
system. It is quite reasonable to assume that one damaging event occurs for each cycle of 
the oscillation. Unfortunately, when the system has multiple frequency components in its 
response, as in a multi-mode structure, the task of estimating the rate of damaging occur- 
rences is very difficult. In the present section an example is given of the classical approach 
to estimate the damage rate in a beam with multiple resonant modes. Improvements on 
this method will be discussed in subsequent sections. 

The main difficulty in estimating fatigue life in multi-mode systems is illustrated in 
figures II. 2.1 through II. 2. 3. These figures show predicted stress in the time domain for a 
beam having one, two, and three resonant modes. The spectra of the stress response for 
each case are shown in figure II.2.4. Figures 11,2.1 through II.2.3 show that as the number 
of modes in the structure is increased, the time domain response becomes increasingly 
complicated. If we were to define a damaging event to be the occurrence of a positive 
stress peak, or where the stress reaches a maximum value in time and then decreases, then 
it is clear from the figures that when th ree modes contribute to the response the rate at 
which damage is accumulated is greater than in the single mode system. 

While it is clear that the time domain response is profoundly influenced by the n umb er 
of modes in a system and that this should in some way affect the fatigue life, it is not 
a simple matter to determine exactly how the damage rate increases with the number of 
modes. When the time history is complicated as in.figureTL5l3, we must make assumptions 
about what causes damage and then attempt to estimate the rate at which the damaging 
event occurs. The problem of selecting a proper definition of a damaging event has not 
received sufficient attention in recent acoustic fatigue research. This will be addressed in 
more detail in the following sections. 

In the current section, applications of two different definitions of damaging event are 
investigated; one in which damage is assumed to occur for each stress peak and one where 
damage results only from the highest positive stress peak between each zero crossing. For 
a linear beam with multi-mode response, approximate analytical solutions for the damage 
rate and fatigue life are developed and compared to estimates obtained using numerical 
simulations for each assumed definition of damaging event. 

The structure studied here consists of a base excited beam with clamped boundaries. 
The ends of the beam are thus assumed to have zero slope and prescribed random trans- 
verse displacement. The axial deflections at the beam ends are assumed to be zero. The 
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assumption of linear response allows us to express the forced vibration in terms of a super- 
position of the eigenfunctions and natural frequencies of the beam. The equation governing 
the transverse deflection, W(x,t ) is 


pAW + El— + rjW= - P AW 0 (t), (II.2.1) 

where p is the mass density, A is the cross section area, I is the moment of inertia, E is 
Young’s modulus, and 77 is a viscous damping coefficient. W 0 (t) is the random acceleration 
that is prescribed at the ends of the beam and is assumed to be Gaussian with zero 
mean and a constant single sided power spectral density G ^ with units of (Wo) 2 /Hertz. 
Equation (II.2.1) is derived in section (II.2.4) and is equivalent to equation (II.4.23) with 
the nonlinear coefficient, c(t) set equal to zero. To solve for the response by modal analysis, 
let 

CO 

W(x,t) = ^2 a i(t)<t>i(x), (II.2.2) 

1=1 

where the <f>i{x) axe the beam eigenfunctions and or, -(f) are unknown functions of time. 

For a beam with clamped ends, the eigenfunctions are given by 


4>i{x) = cos (pix/l) — cosh(p,x//) + Di(sin(pix/l) — sinh (pix/1)) (II.2.3) 


where p, and Z?j are given in the table below for the first six eigenfunctions. 


i 

Pi 

Di 

1 

4.730040745 

-0.982502215 

2 

7.853204624 

-1.000777312 

3 

10.99560784 

-0.999966450 

4 

14.13716549 

-1.000001450 

5 

17.27875966 

- 1.000000000 

6 

20.42035225 

- 1.000000000 


Table II. 2.1 Eigenfunction Coefficients for a Clamped Beam. 

Substitution of equation (II.2.2) into (II.2.1) leads to a set of equations for the un- 
known <*,-(<), 


where 


+ (ai+u;fa t = FiW 0 (t ), i = 1, 2 , . . . , oo, 


2 El fPi \ 4 ,, bh 3 

u ‘ = 7a[j) ’ {I = l2 ) ’ 


c = 


(II.2.3) 


pA' 

1 f 1 

Fi = — y j <f>i{x)dx. 
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b is the beam width and h is the thickness. The solution of equation (II. 2. 3) can be 
obtained as 

F‘ [* 

a «(0 = H / sin(u;rf,.(< — t))W 0 {t)<1t (II.2.4) 

Udi Jo 

where u 2 d . = u> 2 —( 2 / 4, and a°(t) is the transient solution due to a non-zero initial condition. 
Since the transient solution decays in time, it has no effect on the steady-state response. 
We shall then take a°(t) = 0 in the subsequent discussion. The normal strain of the beam 
in x direction is given by 


= _ 

dx 2 




(II.2.5) 


The maximum strain occurs at the top or bottom surface at each end of the beam, y = h/2 
and x = 0 or /. Let 

h e 2 <j>i{ o) 


a, = 


2 dx 2 


The maximum strain is written as 


OO 


(II.2.6) 


(II.2.7) 


The maximum strain rate is given by 


i=l 


OO 




(II. 2. 8) 


t=l 


Equations (II.2.7) and (II.2.8) are often truncated to a finite summation in real computa- 
tion, leading to 

N N 

e = afttiCt), e = «*•«<(*)» (II.2.9) 

•=i i=i 

where N is a finite integer. 

It can be shown that the maximum strain e and strain rate e given by the truncated 
series are stationary Gaussian processes. The first and second order statistics of e and e 
completely determine the process, and can be obtained as follows 

E[e] = E[e] = 0, 

N N 

«=1 j= 1 

N N 

E ^]=EE a ‘“> E [“<( ( )M0], 

i=i i= i 

N N 

eh = EE“.**(«KW], 


(II. 2. 10) 


where 


F F G - f* 

E[aj(<)aj(<)] = — — - — — / e~^ T smuj d .T sinu d rdr, 
2u> di u dj J 0 


E[A, (()»/()] = 


Wfr. /* -< r 


r-{ 


c 




2u d 

mo* '• 


/ '~ <T {si 
2 a>d,. y 0 l 


C 

cos aJd ; r cos u;^ t — cos u di r sin u dj r 

2u > d j 1 

c 2 .1 

- sm u d . r cos r + sin u d . r sm . r > dr, 

. fodiUdj J 


— sin Wd,. r sin u dj r 1 dr, 
2a>d. J 


sm Wd,. r cos r 


and E[-] is the expected value. It can be shown that in the steady-state as t — ► oo, we have 


N N 


<*- E[s 2 ] = EE 




“( [C 2 + (w* - Wrf, ) 2 ] [C 2 + (u>di + Wdj ) 2 ] ’ 


N N 


^ = e[« 2 i = EE 


Ca.a,r i r,-g^(u.? +^S) 


(II.2.11) 


2 [C ! + (wj, - wj, ) 2 ] [C 2 + (w* + ^d, ) 2 ] ’ 


and 

E[ee] = 0. 

Having obtained the statistics of the strain and strain rate, e and e, we then proceed 
to estimate the fatigue life of the system. There may be different approaches to this prob- 
lem dependent on the counting scheme for damaging events that one chooses in applying 
Miner’s rule as discussed in Section II. If we count all the peaks of the strain (or stress) 
as damaging events, and use the fact that e and e are joint Gaussian variables, we can 
take advantage of existing methods as presented in reference [3]. To apply the theory, it 
is necessary to evaluate two parameters: E[Mt], the average rate at which strain peaks 
occur, and a parameter, a, which is the ratio of the rate at which zero crossings occur to 
the rate of strain peaks. The parameter a appears in the probability density function of 
the strain peaks. E [Mr] and a are given by 


where 


E [M t ] = 


<73 

2tt<72 ’ 


a = 


<7102 


/ oo 

OO 

/ OO 

•OO 

fOO 

a\ = / u> 4 $ e (u>)du, 

J — OO 


(II.2.12) 


(II. 2. 13) 
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(II. 2. 14) 


$ £ (u;) is the two sided power spectrum of the strain, and is given by 


N N 
»=1 j= 1 

where 

G ■ F F' 

= ”4^ (II.2.15) 

and 

* i(a,) = ^-J+,cu, ( ' 2 = - 1 ’- < n - 2 - 16 ) 

*v( w ) is the cross power spectrum of a* and aj. In the present study, we assume that 
the random excitation W 0 (t) is bandlimited with a cutoff frequency / c . f c is taken to be 
1000 Hz in the numerical study reported next. This bandlimited assumption is necessary 
for the convergence of the integral for a\ in equation (II.2.13). 

The probability density function for response peaks is given by 


Ppeak{t) 


(1 -a 2 )2 e _ e 2 [2g 2 (1 _ Qi)] -i 
\/2Trai 


ae 





(II. 2. 17) 


where a and o x are defined in equations (II.2.12) and (II.2.13). 

In a random process in which there is only one peak for each zero crossing, a will equal 
to unity and equation (II.2.17) reduces to a Rayleigh density. This is typically considered 
a narrowband process. A process with an infinite number of peaks for each zero crossing 
will have a equal to zero and equation (11.2.17) will reduce to the Gaussian density. The 
process is then said to be broad-band. 

To calculate the damage rate as in equation (II.2.7), m(S) is replaced by a constant 
E [Mt] and equation (II.2.17) is substituted for the peak density. Since the damage rate 
m equation (II.2.7) is expressed in terms of stress, S, rather than strain, e, we must also 
employ the relation S = Ee, where E is Young’s modulus for the material. The damage 
rate in terms of an integration over the strain then becomes 


E (A] - f E [^r]^^-p P eafc(e)de, (II.2.18) 

J — OO C 

with E [Mt] given in equation (II. 2.12) and p pea jt(e) given in equation (II. 2. 17). 

The damage rate as estimated by equation (II. 2. 18) is based on the assumption that 
damage occurs every time the strain reaches a maximum and then decreases. An alter- 
native to assuming that each peak causes damage is to assume that the damage rate is a 
narrowband process. The expected value of the damage rate, E[A] may then be estimated 
by approximating the strain as a narrowband process. The peak density then corresponds 


10 



to the Rayleigh density and the rate at which peaks occur is the same as the rate of zero 
crossings, n 0 (e), 


n 


.(«) = 




2ira f 2 tt(Ti 

The average damage rate and fatigue life are then obtained as in equation (II. 2. 12) 


(II.2.19) 


e Pm] = =£>(VWr (t±i) , 


(II. 2. 20) 


and 


e [D m y 


where as = Ea f , and E is the Young’s modulus of the material. 


(II. 2. 21) 


In the numerical results reported herein, the number of terms used in the expansions 
for the response, equations (II.2.9), is taken to be N = 6. Because of the symmetry of the 
structure, the generalized forces associated with the even modes in equation (II.2.3) are 
zero, 

F 2 =F 4 = F 6 = 0. (II.2.22) 

This indicates that the modes that are anti-symmetric about the middle point (x = 1/2) of 
the beam do not contribute to the long-term dynamic response, and hence to the fatigue 
damage accumulation. 


The methods presented above have been applied to estimate the fatigue life of the 
clamped beam as a function of the level of the random base excitation and as a function 
of the fatigue exponent, 6, in equation (II. 2.1). The predictions were performed using 
the analytical methods discussed above and also by a direct numerical simulation of the 
random response in the time domain. 


Figure II.2.5 shows a comparison of fatigue lives obtained by the analytical method 
discussed above and by numerical simulation. In this case it is assumed that each stress 
peak in the time domain causes damage. The figure shows that the numerical and analytical 
methods give essentially identical results. The effect of including one, two, or three modes 
in the model is also shown. It is found that with the assumption that each stress peak 
causes damage, increasing the number of modes tends to reduce the fatigue fife. It could 
be argued that in general, one obtains a nonconservative fatigue life estimate by only 
accounting for the lowest frequency mode in a structure. 

The calculations shown in figure II.2.5 are repeated in figure II.2.6 except that in this 
case it is assumed that only one damaging event occurs for each time the stress crosses 
zero with positive slope. Here it is assumed that the damage results only from the highest 
peak between zero crossings. It is again found that the analytical and numerical methods 
are in excellent agreement. As in figure II.2.5 increasing the number of modes in the model 
decreases the fatigue life. 

Comparing the results in figures II.2.5 and II.2.6 shows only very small differences 
between the two damage counting schemes. The close agreement between the elementary 
counting schemes examined here should not be interpreted as evidence that all damaging 
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event definitions will give identical results. A more detailed comparison of damage counting 
schemes will be presented in the following sections. 


II.3. Effect of Nonlinearity on Fatigue Life (Peak Counting) 

In the special case where the response of a nonlinear structure is dominated by a single 
resonant mode, the methods discussed above can be extended to estimate the fatigue 
life including nonlinear effects. A single mode model of the displacement response of a 
base-excited beam with large deflection may be expressed as 

a + u 2 0 (a + £a 3 ) + fia = f(t), (II.3.1) 

where a is the response of the first mode so that the beam displacement, W(x , /), is 

W(x,t) = a(t)<f>(x), (II. 3. 2) 

is the corresponding eigenfunction. io a is the natural frequency when the deflections 
are small, /3 is the damping coefficient and e is a constant which determines the degree 
of nonlinearity. f(t) is a Gaussian random process. The derivation of equation (II.3.1) is 
presented in the following section. 

The nonlinearity in the beam response is assumed to result from the stretching of the 
fibers along the beam’s neutral axis which occurs when the deflections are large. When 
large deflections axe included in the analysis the stress (or strain) is not linearly related 
to the displacement. For the nonlinear beam with a single resonant mode, the maximum 
stress (at the beam end) is found to be 

5 = £ (" 2 K + ' ! (() 2 ")’ (II.3.3) 

where E is Young’s modulus, l is the beam length, h is the thickness and d \ , and p are 
constants given by = 6.1513 and p = 4.730. 

For a single mode system described by equations (II.3.1) and (II.3.3), it is possible to 
apply equations (II.1.7) and (II. 1.8) in a fairly rigorous fashion. The result is a nonlinear 
analog of Miles’ linear single degree of freedom analysis shown in equation (II. 1.12). It 
should be emphasized, however, that real structures very rarely behave as single degree 
of freedom systems and that fatigue life estimates are affected by the number of modes 
included in the analysis. The main purpose of the present study is to develop accurate 
and practical methods of applying equation (II. 1.7) to multi-mode nonlinear systems. 

For the nonlinear system considered here the easiest approach is to express the ex- 
pected value operation in equation (II. 1.7) a s an integration over the modal response, o, 
rather than over the stress, S. This leads to 

E [A] = f ra(£(<*)) ^ Q ^ p P eak(at)da. (II.3.4) 

J — oo c 
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It is not a simple matter to derive an exact expression for p pea k(<x ) or m[S(a)) in 
equation (II. 3.4) when the input to the system is Gaussian white noise. If it is assumed 
that there is only one peak for each zero crossing with positive slope (which implies a 
narrow band process) then p pea k ( a ) and m may be calculated so that equation (II.3.4) 
becomes 


E[A] = -w\ (°° |S(a)l*(<* + ia 3 )e 
C J-O O 



(II.3.5) 


where 



(II.3.6) 


and where <r^ o is the mean square response of a(t) when e is set to zero in equation (II.3.1). 

The calculation of the damage rate in equations (II. 3. 5) and (II. 3. 6) relies on knowledge 
of the peak probability density which for the present system is given by 


Ppeak(&) 


(a + ea 3 ) 



(II.3.7) 


This expression may be obtained from the joint probability density of a and ct if 
the response is a narrowband process. Unfortunately, the joint density of a and d may 
be obtained only for rather specialized nonlinear systems. Since our goal is to develop a 
method that is applicable to a wide class of systems, an approximate method of estimating 
the damage rate will now be described that includes the effects of nonlinearities on a 
through the method of equivalent linearization [4]. 

The method of equivalent linearization consists of replacing equation (II.3.1) with an 
equivalent linear system given by 

a + ula + /3d + e(a) = f(t), (II.3.8) 

where u 2 t is the equivalent linear frequency and e(a) is an error term. is chosen so 
that the mean square of the error is minimized. If the error, e(a), is then neglected, 
equation (II. 3. 8) is a simple linear system with a natural frequency that depends on the 
spectrum level of the excitation. Since the system is linear, its response will be Gaussian 
for a Gaussian input, f(t). This leads to 

<■>; = ^(1 + ^1 + 1210. (II.3.9) 


then 


If the single sided power spectral density of f(t) is G /(£-) with units of f 2 / Hertz, 


_ G f 

^Ot 0 


4.-r < IL31 °) 

where Gf is assumed to be a constant at all frequencies. One may then calculate the mean 
square response of the equivalent linear system to be 


2 g / 

“• 4U&0* 


(II.3.11) 
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With the response of the equivalent linear system assumed to be Gaussian when e(a) 
is neglected in equation (II.3.8), one may calculate the peak probability density to be a 
Rayleigh distribution as in equation (II. 1.10). The peak density for the equivalent linear 
system is then 

a — 9 4- 

PM = — e , (II. 3. 12) 

°a c . 

where cr£ e is given in equation (II.3.11). 

To evaluate equation (II.3.4), we must also approximate the rate of occurrence of 
peaks, m(5(a)). For the present study this will be taken to be a constant equal to the 
equivalent linear natural frequency of the oscillator in Hertz, u} e /2n. 

By using equivalent linearization, equation (II.3.4) may then be approximated by 


E 1 A ' = I 





(II.3.13) 


where 5(a) is given in equation (II.3.3). 

The comparison of equations (II.3.5) and (II.3.13) may be simplified somewhat by 
defining a normalized random variable, 


7 = 



(II. 3. 14) 


Substitution of equation (II. 3. 14) into equations (II.3.5) and (II.3.6) gives 
E[A]* = U ° Jo°° l‘ S (^ cr °)i b ' (7 + d^y 




C 


.-il-irr 2 




(II.3.15) 


where the R subscript on E[A] denotes the ‘rigorous’ solution. 

Substitution of equation (II.3.14) and the use of equations (II.3.9) and (II.3.10) in 
equation (II.3.13) give 


E ^ E = ^J 0 I‘ S '(7^o)| 6 7(^(1 + \/T+ 12e<72)) 3/2 e -^( 1+ \/ l + 12 ^. 2 )d 7 , (II. 3. 16) 

where the E subscript on E[A] denotes the equivalent linearization solution. 

Substituting equation (II.3.14) into (II. 3 . 3) and rearranging give 

•S' = aia 0 (7 + C7 2 ), (II.3.17) 

where 

ai = E/i(j) 2 , (II.3.18) 
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and 


<r 0 di 



(II. 3. 19) 


The results of this section may be used to estimate the effects on fatigue life of nonlin- 
ear strains in a single mode system and to assess the applicability of equivalent linearization 
to the nonlinear fatigue problem. Figures II.3.1 and II.3.2 show predicted fatigue lives of a 
base excited beam for a range of excitation levels and for the damage exponent, 6, in equa- 
tions (II.3.15) and (II. 3. 16) equal to 4 and 6. The geometrical and material parameters 
of the beam are identical to those used in the results presented in Section II.2. Excellent 
agreement is shown between the results of equations (II.3.15) and (II.3.16) for both values 
of b. 

Figures II.3.1 and II.3.2 show the predicted fatigue life with and without the effect 
of nonlinear axial strain. In the calculations, the nonlinear axial strain was neglected by 
setting a = 0 in equation (II.3.17). In this case the nonlinearity affects the displacement 
response through equation (II.3.1) but the relation between strain and displacement is 
assumed to be linear. This is essentially the same approach taken in reference [4]. The 
figures show that when the full nonlinear strain-displacement relation is used, the added 
stretching strain substantially reduces the fatigue life at high excitation levels. 

The results shown in figures II.3.1 and II.3.2 indicate that nonlinear stiffness in the 
beam can reduce the response amplitude relative to that of a linear system and hence 
increase fatigue life. However, the figures also show that the nonlinear stretching strain 
results in an increase in the damage rate and hence partially cancels the beneficial effect 
of the reduced displacement. 


11*4 Governing Equations For A Nonlinear Beam 

The governing equations for a beam will now be developed including what are likely to 
be the dominant nonlinear effects. The nonlinearities accounted for here are due to axial 
stretching and friction damping at the boundaries. The ends of the beam are assumed to 
be clamped in a fixture that allows axial slipping when the axial force exceeds a certain 
value. For simplicity, this slipping is assumed to occur at only one end, the other end 
being held rigidly. 

The friction damping model has been included in this analysis because it is suspected 
to be the most likely nonlinear damping mechanism in the beam. When the experimental 
phase of the current project is complete, data will be available to determine the relative 
influence of stiffness nonlinearity (due to axial stretching) and damping nonlinearity (due 
to friction at the boundary). Although the friction damping model is developed in this 
section, it has not been utilized in a numerical model to assess its influence on fatigue life. 

The derivation of the equations of motion will be based on Hamilton’s principle, 



- V + W)dt = 0, 


(II.4.1) 
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Where S is the variational operator, T is the kinetic energy, V is the potential energy and 
W is the work done by nonconservative forces. The limits of integration, ti and t 2 define 
an arbitrary time interval. 

If the beam is subjected to a base acceleration, W 0 (i), and if the deflection of the 
beam mid-line relative to its moving supports is W(x,t), where x is the axial coordinate, 
then the kinetic energy may be written as 

T=^-J\w 0 (t ) + W(x,t)) 2d x, (11.4.2) 

where p is the density, A is the constant cross sectional area, l is the beam length, and the 
dots denote partial differentiation with respect to time, t. It is assumed that the beam is 
thin so that rotary inertia effects are negligible. 

The potential energy will be assumed to consist of contributions due to strain energy 
in the beam, V s , and the energy V&, stored in a linear, axial spring at the boundary at 
x = l. To construct the strain energy, V a , it will be assumed (as in elementary beam 
theory) that the strain is dominated by the axial strain, e ZI . This leads to 


V, 


=UL 


Et\ x dzdx , 


(II. 4. 3) 


where b is the beam width (assumed to be constant), h is the thickness, E is Young’s 
modulus, and z is the transverse distance relative to be beam mid-line. 

Bringing the variational operation inside the integrand of equation (II.4.1) gives 


f * (6T-6V + 8W)dt = 0, 

Jt i 

where from equation (II.4.2), 

ST = pA [ (W 0 (t) + W(x,t))8W(x,t)dx, 
Jo 


(II.4.4) 


(II. 4. 5) 


and, from equation (II.4.3), 

6V = SV 3 + SV b 




Ee xx 8e xx dzdx -|- 


(II. 4. 6) 


In equation (II. 4. 5) we have used the fact that W 0 (t} is prescribed so that its variation is 
zero. The contributions due to nonconservative forces in SW will be discussed later. 

Substituting equation (II. 4. 5) into equation (II.4.4) and carrying out the integration 
of ST over time by parts give 



ST dt = —pA 



(W 0 (t) + W(x,t))6W(x,t)dxdt, 


(11,4.7) 
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where the variation of W(x, t ) is zero at t = t \ , and t = t 2 . 

The variation of the strain energy in equation (II.4.6) must be written in terms of 
convenient displacement coordinates. This is accomplished by expressing the axial strain 
in terms of U(x,t), which is the axial displacement at 2 = 0, and the transverse 

displacement at 2 — 0. This gives 


_ dUjxJ) d^Wjx^t) 1 f dW(x,t) \ 2 

dx x 2 2 \ dx ) 


(II. 4. 8) 


The last term on the right hand side of equation (II.4.8) accounts for axial stretching which 
will be important when the deflection is large. 

Substituting equation (II. 4. 8) into equation (II.4.6) gives 


6V, = bE 



- ^rr + 2 Wl){6U x - z6W xx + W x 6W x )dzdx, 


(II.4.9) 


where the subscripts denote partial differentiation. Integrating equation (II. 4. 9) over z 
gives 

6V, = bE J\^W XX 6W XX + h 2 (U x W x + \wl)8W x {U x + \w 2 x )dx. (II.4.10) 

The finite stiffness of the beam support in the axial direction may be modeled by an axial 
spring at x = l having spring constant, k. The potential energy of this spring is 

H = \kU 2 {l). (IL4.11) 

Integrating equation (II.4.10) by parts and including the variation of the boundary poten- 
tial energy in equation (II.4.11) lead to 


8V = Eb 


h 3 
12 


W T 


- + \wl)W^6Wdx - + \wDiVdx\ 

+Eb {£(W„6W x |J - W IX ,6WI' 0 ) + h(U,W, + 


+Ebh(U x + hvl)SU\[ + kU(l)6U(l). 

(II.4.12) 

The remaining term to be evaluated in Hamilton’s principle in equation (II. 4. 4) is the 
contribution due to nonconservative forces. This includes the effects of applied forces and 
damping. It will be assumed that there are no external forces acting other than those that 
produce the base excitation which is already accounted for in the kinetic energy expression, 
equation (II.4.2). Two different damping mechanisms will be included in the model. One 
will be viscous damping which depends only on the velocity at each point along the beam 
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length and the other is due to friction at the support at x = l. The model of frictional 
damping accounts for sticking and slipping depending on the axial force at the boundary. 
This results in hysteresis which greatly complicates the analysis of the random response. 

The variation of the work done by viscous damping (the virtual work) may be modeled 
as 

= - j rfW(x,t)8W{x,t)dx , (II.4.13) 

Jo 

where r] is the viscous damping coefficient. The negative sign accounts for the fact that 
damping consists of the removal of energy from the system. 

Frictional damping at the support will be modeled as a contact surface with a clamping 
force N. It is assumed that as long as the magnitude of the axial force is less than a 
critical value, fiN , there is no slipping and the axial force applied by the beam is equal and 
opposite that of the frictional support. As the axial force at the beam end is increased to 
fiN, however, the contact surface slips and no further increase in axial force is allowed. The 
stick/ slip boundary may then be thought of as a force limiter. The behavior is characterized 
in three possible states, (1) axial force equal to -fiN for slipping in negative x direction; 
(2) zero velocity (no slip) when the magnitude of the axial force is less than fiN, and (3) 
axial force equal to +fiN for slipping in the positive x direction. 

To construct the equations that describe the axial boundary condition let the frictional 
damping force be denoted by f s (U(l),U(l)). The virtual work associated with friction 
damping is then 

SWf = f a (U(l), U(l))8U(l). (II. 4. 14) 

The axial boundary condition at x = 1 follows from letting 8W = 8W V + 8Wf, then 
substituting equations (II. 4. 12) through (II.4.14) into equation (II. 4. 4) and combining the 
terms that multiply 8U(1). Setting the result to zero, we have 


-Ebh(U, + 


kV(l) + f.(U(l),V (‘ )) 


su(i) = o. 


(II.4.15) 


The first term in equation (II. 4. 15) may be evaluated using the governing equation for 
axial deflection which may be obtained from equations (II.4.2), (II.4.4), and (II. 4. 12), 

Ebh-^(V, + = 0. (11.4.16) 

This equation contains no inertia terms because axial kinetic energy was neglected in 
equation (II.4.2). It is assumed that resonances in the axial direction are outside the 
frequency range of interest. 


Integrating equation (II. 1.16) over x gives 


U x + -Wl = c(t), (II. 4. 17) 

where c(<) depends only on time. From equations (II.4.8) and (II.4.17) the axial strain is 

d 2 W(x,t) 


£ n — c(t) — z - 


dx 2 


(11-4.18) 
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where the first term c(t) is seen to account for axial stretching of the raid-line at z — 0, 
and the second term is the usual bending strain. The stretching of the beam mid-line, c(t) 
is a constant along the beam length. 

Substituting equation (II.4.17) into equation (II.4.15) leads to 

-Ebhc{t) - kU{l) + *7(0) = 0. (II.4.19) 

It is convenient to express this boundary condition in terms of the axial stretching, c(t), 
and the transverse deflection, W(x,t), by eliminating * 7 ( 0 - This may be accomplished by 
integrating equation (II. 4. 17) over x, 

f l 1 

U ( x > 0 + J o = xc ( f ) + «(*). (II. 4. 20) 

where ( is a dummy variable and a(t ) is a constant of integration. Since the axial dis- 
placement is constrained to be zero at x = 0, i.e., U(0, t) = 0 then from equation (11.4.20), 
a(<) = 0. Evaluating equation (II.4.20) at x = / then gives 

V(l,t) = lc(t)-lj'wl((,t)d(. (II.4.21) 

Substituting equation (II.4.21) into equation (II.4.19) and rearranging give the friction 
force in terms of c(t ) and W(x,t) 

f s (c(t), W(x, 0) = (kl + Ebh)c(t ) - | j 0^- (II.4.22) 

Equation (II.4.22) must be solved along with the governing equation for Il r (x, t) which 
follows from equations (II.4.4), (II.4.7), (II.4.12), (II.4.13) and (II.4.17), 

£bh z 

pAW(x,t ) + — -W xxxx - Ebhc(t)W xx +t]W = - P AW 0 . (II.4.23) 

Since IT(x, t) is the transverse deflection of the beam mid-line relative to the motion of 
the supporting clamps, it will be assumed that W(x,t) and H7,.(x,<) are zero at x = 0 and 
x = l. This eliminates the boundary terms for W(x,t) in equation (II.4.12). 

Equations (II.4.22) and (II.4.23) must be solved simultaneously for c(t) and W(x, t) 
with the friction force f 3 properly accounted for. Due to the nature of the friction force, an 
iterative solution procedure must be used. In the present study, the response is simulated 
in the time domain using Newmark’s numerical integration scheme. 
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Figure II.2.4 Predicted stress spectra of a linear beam. 
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Figure II. 2. 5 

Predicted fatigue life of a linear beam with multiple modes. 
Damage is assumed to occur for each stress maximum. 
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Figure II. 3.1 Predicted fatigue life of a nonlinear beam with one resonant mode. 






III. 1 Rainflow Counting 


The methods used in the previous sections for identifying damaging events are appeal- 
ing because they require statistics of the random response that may be easily computed 
using exact closed form expressions. In the present section we estimate damage based on 
a consideration of the mechanical behavior of materials rather than on easily computed 
statistics. While the resulting procedure does not permit simple analytical expressions to 
be used to estimate fatigue life, it has been shown to produce more accurate estimates of 
fatigue life than the methods discussed above. 

A great deal of research on fatigue is concerned with estimating the usable lifetime 
of structures that are subjected to complicated loading. These studies are usually broken 
into two cases, one where the structure experiences significant plastic deformation and 
one where the deformation is predominantly elastic. When the strains are sufficiently 
large to induce plastic deformation, the damage mechanisms are different from the case 
where elastic behavior predominates. Since the structure lasts considerably longer when 
the strains are elastic, this situation is referred to as high-cycle fatigue. In studies of 
acoustic fatigue, the strains are considerably smaller still and damage accumulates very 
slowly. Acoustic fatigue is a special case of high cycle fatigue. 

In contrast to plastic, or low-cycle fatigue, in which damage can be directly measured 
by observing the growth in crack length, the damage in high-cycle fatigue is microscopic 
and during the majority of the structure’s life, it can only be measured using sophisti- 
cated techniques. It is assumed that unobserved damage accumulates over the life of the 
structure. The time duration for this damage to result in failure is the fatigue life. 

Although the strains in high-cycle fatigue Eire primarily elastic, damage is always 
associated with inelastic behavior. It must be assumed that while the majority of the 
material responds elastically, there must be high strains at specific locations within the 
structure that result in microscopic damage. In acoustic fatigue studies, the deflections, 
stresses, and strains are estimated with the assumption of elastic response but the damage 
estimation should be based on considering the strains to be plastic. 

It is well known that when a material experiences cyclic plastic deformation, a hys- 
teresis loop is observed in a plot of stress versus strain. While it may not be possible 
to measure the stress and strain at microscopic areas where damage occurs in high-cycle 
fatigue, it is reasonable to assume that a hysteresis loop exists. This is the basis of modem 
damage counting schemes such as Rainflow cycle counting. 

Rainflow cycle counting was first proposed by Matsuishi and Endo in 1968 [1], Prior to 
this development considerable effort had been made to construct a cycle counting scheme 
that correlated with a physical understanding of the damage process and with experimental 
results. 

The most often cited work on acoustic fatigue prediction is that of J.W. Miles [2] 
as discussed in section II. In this paper a remarkably simple formula is presented for 
estimating the fatigue life of a simple panel subjected to random loading. It is assumed 
that the strain response is dominated by that of the first resonant mode so that it may 
be considered to be a narrowband process. This narrowband assumption greatly simplifies 
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the estimation of fatigue life because it may then be claimed that whatever the damage 
mechanism, there will be a damaging event for each cycle of the oscillating strain or stress. 
By assuming that the rate at which damage occurs is equal to the oscillation frequency, 
one can use the probability density for stress maxima to compute the expected value of 
the fatigue life. The details of the calculation are given in references [2-4]. 

Although the response of complex structures can rarely be considered a narrowband 
process, there has been considerable motivation to assume that damage is easily related 
to some statistic that can be readily calculated. If one assumes that damage occurs each 
time the stress or strain experiences a maximum or minimum it is possible to employ a 
similar formulation to that of Miles [2] to estimate fatigue life. An expression for the 
probability density of strain peaks in a wide-band random process is given in [4]. This 
result depends on the ratio of the expected number of zero crossings at positive slopes to 
the expected number of peaks. This quantity is easily computed from the response power 
spectral density. An example of the application of this approach is given in the previous 
sections and is also presented in Vaicaitis and Choi [5]. 

While the assumption that damage occurs for each strain peak has been primarily 
motivated by its mathematical utility, there have been several attempts to justify it ex- 
perimentally. Broch [6] conducted an extensive fatigue test using a total of two hundred 
samples. Half of these samples were driven such that response behaved like a one degree of 
freedom system and the remaining samples responded as two degree of freedom systems. 
An attempt was made to control the experiments so that the mean square response levels 
and the total number of zero crossings were the same for each sample. In the one degree 
of freedom tests the ratio of the number of zero crossings to the number of peaks, a, was 
approximately unity and in the two degree of freedom tests the ratio was 0.275. An at- 
tempt was made to correlate the data using the known statistics of strain or stress peaks. 
It was concluded that further work was necessary before firm conclusions could be drawn. 

Hillberry [7] conducted a similar comparison of the effect of response bandwidth where 
one set of data was obtained using narrowband response and the remaining data was 
broadband. In his broadband data , however, the ratio of the number of zero crossings to 
the number of peaks was 0.79 which is not significantly different than unity as in the case 
of narrowband data. Only minor differences were detected in the fatigue lives due to the 
different loadings. In a later work, Linsley and Hillberry [8] conducted fatigue tests with 
loadings that had significantly different spectra but identical peak probability densities. 
It was found that the fatigue lives were very similar even though the response spectra 
were quite different. They concluded that the peak probability density contained all the 
necessary information for fatigue predictions. 

While Linsley and Hillberry’s result [8] is interesting and lends support to the claim 
that peak counting may be sufficient for fatigue prediction, it is found that two signals that 
have the same peak probability density could also have similar damage rates as predicted 
using Rainflow counting. This was observed by Wirsching and Shehata [9] who estimated 
fatigue lives using Rainflow counting on simulated data sets that had eleven different power 
spectra. These spectra had a ratios that varied from 0.218 to 0.998 . It was found that 
the mean and standard deviations of stress cycles as determined by Rainflow counting are 
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directly related to a. The authors conclude that for the purpose of estimating fatigue life, 
the spectral shape of the response can be characterized by the total power and a. Further 
justification of this claim is needed. 

Since the early seventies, there has been general agreement that for complicated ran- 
dom loading it is necessary to use a scheme such as Rainflow counting to obtain reliable 
fatigue estimates. This type of procedure is especially helpful in properly accounting for 
the effects of mean stress superimposed on the fluctuating loads. The mean stress could 
be the result of very low frequency oscillations such those caused by thermal stress, geo- 
metrically nonlinear in-plane strains due to large deflections, or low frequency structural 
modes. The effect of mean stress could be very important in acoustic fatigue analysis. 

The manner in which Rainflow counting accounts for mean stress is illustrated in figure 
III. 1.1 which is taken from reference [10]. The key to the method is to attempt to count 
closed hysteresis loops. As discussed above, although these stress-strain hysteresis loops 
may be unobservable due to the difficulty of measuring stress and strain at the locations 
of damage accumulation, there must be some hysteresis in the structure because damage 
will never occur if the response is perfectly elastic. Hysteresis loops associated with two 
rather similar strain time histories are shown in figure III. 1.1. In each time history a 
single large amplitude strain cycle occurs before the strain oscillates at a lower amplitude. 
The difference in the two time histories shown is that the half cycle preceding the smaller 
amplitude oscillation is negative in history A and positive in history B. As indicated in the 
figure, the hysteresis loops corresponding to these cases show that for history A, during 
the small oscillation cycles there is a positive mean stress, er 0 , and for history B the mean 
stress is negative during the small amplitude cycles. The damage associated with these two 
histories is expected to differ substantially because it is known that a positive mean stress 
can shorten fatigue life and a negative mean stress can lengthen it. It is important to note 
that older cycle counting schemes (such as simple peak counting) that do not consider the 
hysteresis loops will predict identical damage rates for these two time histories. 

To understand the way Rainflow counting accounts for hysteresis we need to describe 
the method in detail. There are many excellent descriptions of the approach that are 
available in the literature. The following is adapted from reference [10]. 

The Rainflow counting procedure is best visualized by viewing the strain time history 
with the time axis pointing downward as in figure III.1.2 (adapted from reference [10]). 
The method gets its name from constructing cycles by imagining rain flowing downward 
with the strain curve acting like a series of roofs. The first step is to rearrange the time 
history so that it begins and ends at the strain value of greatest magnitude in the block of 
data. This rearrangement of the data is justifiable because the response is considered to 
be a stationary random process. A shift in the time axis will not affect the statistics. It is 
assumed that rain begins to flow at the beginning of the trace at point A in figure III.1.2 
and at each strain reversal in the time history. To identify the closed hysteresis loops a 
number of rules are applied to the falling rain. The rain is assumed to continue flowing 
downward unless: 

a. The rain began at a local maximum (peak) and falls opposite a local maximum greater 

than that from which it came. 
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b. The rain began at a local minimum point (valley) and falls opposite a local minimum 
point greater (in magnitude) than that from which it came. 

c. It encounters a previous rainflow. 

Figure III. 1.2 illustrates the application of the method. The dots indicate the initiation 
of rainflow at each reversal in the time history. 

A. Rain flows from point A over points B and D and continues to the end of the history 
(which is point A again) since none of the conditions for stopping rainflow are satisfied. 

B. Rain flows from point B over point C and stops opposite point D, since both B and 
D are local maximums and the magnitude of D is greater than B (rule a above). 

C. Rain flows from point C and must stop upon meeting the rainflow from point A (rule 

c). 

D. Rain flows from point D over points E and G and continues to the end of the history 
since none of the conditions for stopping rainflow are satisfied. 

E. Rain flows from point E over point F and stops opposite point G, since both E and 
G are local minimums and the magnitude of G is greater than E (rule b). 

F. Rain flows from point F and must stop upon meeting the flow from point D (rule c). 

G. Rain flows from point G over point H and stops opposite point A, since both G and 
A are local minimums and the magnitude of A is greater than G (rule b). 

H. Ram flows from point H and must stop upon meeting the rainflow from point D (rule 
c). 

The events identified above may now be combined to identify complete cycles, or 
closed hysteresis loops. The rainflow from A to D and from D to A forms a complete cycle. 
Another cycle is formed by the flow from B to C and from C to the strain level at B (in 
the line connecting C and D). There are also cycles formed between strain ranges E-F and 
G-H. 


The rainflow process identifies cycles that are associated with closed hysteresis loops 
as indicated in figure III. 1.3. This figure shows the loops that correspond to each cycle. 
The strain maxima are each located at a tip of a loop. After the strain goes from A to B 
it decreases to C. As the strain increases from C it reaches the level of B again while on its 
way to D. The Rainflow counting procedure correctly indicated that there is a hysteresis 
loop between B and C which is a complete cycle contained within the larger loop between 
A and D. 


Once the cycles and associated strain ranges are identified the, Palmgren-Miner linear 
damage accumulation rule may be applied [11], This rule states that we may simply add 
the incremental damage (or life used up) due to each cycle to estimate the fatigue life. The 
damage due to one cycle may be expressed as 



1/2AS* 

1 - s mi /s u 


(III.1.1) 


where AS, is the stress range of the cycle, S m , is the mean stress of the cycle, S u is the 
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ultimate strength of the material, and b and c are constants obtained from a S-N diagram 
for the material with constant amplitude loading. The manner of incorporating the mean 
stress in equation (III. 1.1) is after Morrow [12]. The total damage is simply 

(III. 1.2) 

% 


Failure is assumed to occur when D = 1. 


III. 2 Experimental Verification of Rainfiow Counting 


There have been a number of experimental studies to verify the accuracy of fatigue 
prediction schemes. In the following, we will describe two reports that address issues 
that are relevant to acoustic fatigue studies. The two most commonly cited works on the 
experimental verification of Rainfiow counting axe those of Dowling [13] and a program 
conducted by the SAE Cumulative Fatigue Damage Division [14]. In Dowling’s study 
a total of 83 specimens were tested to failure. His main purpose was to construct a 
fatigue prediction procedure that could account for the effects of the sequence of the strain 
cycles and could also account for fluctuating mean stresses. A number of complicated 
load histories were used. Data is presented which indicates that Rainfiow counting does 
a very good job of handling time histories where small, high frequency oscillations are 
superimposed on top of large amplitude, low frequency oscillations. 

Signals that have small fluctuations added to large low frequency oscillations are 
commonly seen in the response of complex structures where large amplitude low frequency 
modes coexist with an uncountable number of high frequency resonances. Simplified time 
histories that have this feature are shown in figure III.2.1 (taken from Dowling [13]). 
The traces a and b shown in the figure would have identical fatigue lives if simple peak 
counting were used. Dowling showed that the actual fatigue lives of specimens will differ 
substantially with these two types of response and that the difference in fatigue life is 
accurately estimated only if a procedure such as Rainfiow counting is used. 

Tests were conducted by adding two sine waves to obtain stress and strain histories 
as shown in figure III.2.2 a ( again taken from reference [13]). In these tests the amplitude 
of the small, high frequency fluctuation was held constant and the amplitude of the low 
frequency component, Aej, was varied. Figure III. 2.2 b shows that his prediction procedure 
using Rainfiow counting accounts for changes in Aei remarkably well. 

The data of figure III.2.2 correspond to a low cycle fatigue test where the material 
experienced significant plastic deformation. This is apparent in the figure because the 
traces for strain and stress differ substantially. Another set of specimens were tested 
where the strains were predominantly elastic so that high cycle fatigue mechanisms would 
be encountered. Typical stress time histories for these tests are shown in figure III.2.3. In 
this case a high frequency sine wave was superimposed on a lower frequency triangle wave. 
The amplitudes of the high frequency sine wave and the low frequency triangle wave were 
held constant. The fatigue lives were measured and predicted for a range of values of the 
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frequency of the triangle wave. Comparisons of the measured and predicted results are 
shown in the figure. Note that the vertical scale spans only one decade so that the error 
is roughly a factor of two. This is quite good agreement for a fatigue prediction. 

A similar test to that of figure III. 2. 2 was conducted where the maximum stress range 
was kept constant, A<ti = 100 , but the relative amplitudes of the low frequency and 
high frequency components were varied. It was found that as the amplitude of the high 
frequency component was increased, the fatigue life was reduced dramatically. The results 
are shown in figure III.2.4. It should be noted that simple peak counting would not predict 
that the amplitude of the high frequency component would have a major effect. This is 
because the rate at which peaks occur is held constant and the peak amplitude distribution 
is not changed substantially. Rainflow counting, however, does a remarkably good job of 
accounting for this. 

In the results shown in figures III.2.1 through III.2.4 the stress and strain time histo- 
ries essentially contain subcycles with varying lower frequency (mean) fluctuations. It is 
well known that compressive mean stresses reduce damage and that tensile mean stresses 
increase the rate of damage accumulation. Rainflow counting enables one to account for 
these two effects for each cycle identified by the procedure as indicated in equation (III. 1.1). 

A number of other experiments were reported by Dowling [13] which will not be dis- 
cussed here. These tests included other complicated time histories. He did not, however, 
conduct a comprehensive evaluation of Rainflow counting for very high cycle fatigue prob- 
lems with loadings that axe as complicated as encountered in acoustic fatigue. One may 
conclude, however, that Rainflow counting should lead to vastly more accurate estimates 
of fatigue life than methods that are often used in acoustic fatigue studies. 

Another study that is often cited in support of the use of Rainflow counting is a 
joint effort by the SAE Cumulative Fatigue Damage Division [14]. In this, study three 
different complex loadings were applied to notched steel specimens. A total of 57 specimens 
were tested. The load time histories used are shown in figure III.2.5. The time history 
corresponding to bracket vibration shown in the figure most closely matches that seen in 
aerospace structures. The response power spectra for these histories were, unfortunately, 
not presented in the report. A comparison of predicted and measured fatigue lives is shown 
in figure III. 2. 6. In these predictions the Rainflow counting method was used to identify 
cycles as discussed in Section II. 

An investigation is described in reference [14] of the effect of mean stress on fatigue 
life. Predictions were made with and without mean stress included. This was accomplished 
by either including or neglecting S m in equation (III.l.l). As discussed above, tensile mean 
stress is known to reduce fatigue life while compressive mean stress increases it. It was 
found that for a time history where the average stress was near zero, the reductions and 
increases in fatigue life due to fluctuating mean stresses tend to cancel out. It is likely 
then, that for Gaussian time histories with zero mean, the mean stress effect may be safely 
ignored. When the structure responds nonlinearly, however, such that in-plane stretching 
strains are significant, the mean stress effect could be quite substantial. 

It should be pointed out that the tests conducted by Dowling [13] were conducted 
on smooth samples with uniaxial loading and the tests discussed in reference [14] were 
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conducted on notched specimens in tension and compression. None of the studies reviewed 
to date examined fatigue where the primary mode of deflection is in bending as it is in 
acoustic fatigue. It is possible that this could affect the applicability of these studies to 
the acoustic fatigue problem. 

A large number of other studies have compared experimented fatigue measurements 
with predictions. We have attempted to limit our attention to those studies that relate to 
the acoustic fatigue problem. This means that those studies that investigated low cycle 
fatigue are excluded. Low cycle fatigue is important if acoustic fatigue tests are accelerated 
by increasing response levels to the point where plastic deformation is significant. This 
would drastically change the results and care should be taken when high level tests are 
designed. 
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IV. Power Spectral Density of Multi-Mode Systems 


In the following sections a comparison is presented of fatigue lives predicted using 
simple peak counting to those obtained using Rainflow counting. The results presented 
indicate that when a signal contains significant low frequency oscillations there can be 
considerable difference between the fatigue lives predicted by the two methods. Further 
work is needed to compare the methods in typical aerospace applications. 

The comparisons of fatigue lives estimated using peak counting and Rainflow count- 
ing to be presented in the following were performed using a time domain simulation of 
the response. An attempt has been made to simulate the response of typical aerospace 
structures. The method of simulating the time domain signal is presented in section V. 
The method depends on the response power spectral density. It is possible to accurately 
simulate the time domain response corresponding to any given power spectrum. In the 
present study, power spectra are created to model ‘typical’ structures that have nine reso- 
nant modes. The calculation of the response spectrum of a multi-mode system is presented 
in the present section. By processing the resulting time histories to predict fatigue lives, 
it is then possible to examine the influence of changes in the response spectrum (such as 
modifying the relative contribution of various modes) on the fatigue estimated by the two 
methods. 

In the present study, fatigue lives were calculated for systems having nine resonant 
modes over the frequency range from 0 to 2000 Hertz. As will be discussed in section 
VI, it is found that by increasing the relative contribution of the lowest frequency mode, 
substantial differences are observed between the fatigue lives obtained using peak counting 
and those estimated using Rainflow counting. In every case, peak counting produced 
shorter lives than did Rainflow counting. 

It should be noted that the fatigue estimates using peak counting performed here are 
not based on the single degree of freedom assumption of reference [2] of section III. Fatigue 
estimates using the single degree of freedom approach differed from those obtained using 
Rainflow counting by factors ranging from 48.69 to .0098. This method can be either highly 
conservative or very nonconservative depending on the system’s spectral characteristics. 

The primary result of the present comparison is that Rainflow counting produces 
significantly different fatigue life predictions than simple peak counting. This can also 
be seen by examining the Rainflow counting method. When low frequency oscillations 
are present in a signal along with high frequency components, peak counting will produce 
substantially shorter fatigue lives than Rainflow counting. While this is well known for 
certain time histories, it has not been demonstrated that Rainflow counting is necessary 
for time histories that are typically encountered in acoustic fatigue applications. The 
preliminary results presented here indicate that Rainflow counting is necessary in acoustic 
fatigue studies. 

In the following, expressions are developed for the power spectral density of the re- 
sponse of a system that is subjected to Gaussian white noise excitation and that has a 
number of resonant modes. This will allow one to construct the power spectrum of the 
stress or strain response given the system resonant frequencies, modal damping ratios, and 
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the relative contributions of the modes. The results of this section will be used in section 
V to create power spectra of simple multi-mode systems. 

Suppose the deflection at some point can be expressed as 


W(x,y,t) = (IV.l) 

1=1 


where <f>i(x,y) are eigenfunctions and oti(t) axe modal responses. In a linear system a,(<) 
are solutions to 


Qt -j- 2u?t£fdf| + ol{ — Qii^t ) ? i — 1, 2, . .., oo 


(TV .2) 


We will assume that <7, is some constant, f a (t ) is Gaussian white noise, w, is the natural 
frequency of the mode, and £, is the damping ratio. 

The strain at some point on the structure may be expressed in the form 


OO 

e(t) = ^o l a,(<), (IV. 3) 

1=1 

where the a, indicate the contributions of the various modes. 

The solution for <**(<) in equation (IV.2) is 

t 

= J hi(T r ) qi f 0 (t - T')dT' (IV.4) 

0 


where 

M r> ) = ^ g sinfay/l - £,V), r' > 0. (IV.5) 

It is assumed that oti(t ) = 0 for t < 0. We may compute the power spectral density of the 
strain, e(f), by first computing the strain auto correlation function, 

R (e (r) = E[e(*)e(t + r)], (IV.6) 

where E[-] denotes the expected value. The strain is assumed to be a weakly station- 
ary random process. The power spectral density of e(<) is then obtained by the Fourier 
transform of equation (IV.6), 


*«(w) = 



(IV. 7) 
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Substituting equation (IV. 3) into equation (IV.6) gives 
Ree(r) = E 


r co oo 


a i a j<Xi(i)<Xj(i + T ) 

.=1 j = 1 


(IV.8) 
(IV.9) 
(IV. 10) 

R a{aj (r) = E [ai(t)aj(t + r)] (IV. 11) 

is the cross correlation function for a,(t) and aj(t). Our assumption that the process is 
weakly stationary causes this cross correlation function to be independent of t. Substituting 
equation (IV. 10) into equation (IV.7) and rearranging gives 


= °«' a i E l Qi (0«i(< + r )] » 

*'= i i= i 

OO OO 

= aia,jR Q;ai (r), 

«'= i i= i 


where 




- OO OO °? 

.h = jr E £ *<*> / 
<=1 J=1 -L 


a * a j^a<Oj ( W )’ 

>=1 j=l 

where the cross power spectral density of or,(f) and a,(t) is 

OO 

/ ^ojiry-^dr 

— OO 

Substituting equation (IV.4) into equation (IV. 11) gives 

t t+r’ 

= E (y A.(r') 4i /„((-r')<iT')( J fci(r')®/.(< + r - r')*-') 
0 0 
t 

= 9i«> J J fti(T’)fcy(r")E[/ 0 (i - t')/„(< + t - r")]*-V 


(IV. 12) 
(IV. 13) 


(IV. 14) 


(IV. 15) 


o o 


The excitation, f 0 (t), is assumed to be stationary Gaussian white noise so that 

E[/ 0 (t - T')f 0 (t + t - r")] = 27 r$/ o <S(r + t' - r"), (IV.16) 

where <£(-) is the Dirac delta function and $ y c is the two-sided power spectral density of 
/„(<). Substituting equation (IV. 16) into equation (IV.15) and integrating over t" gives 


T 

RaiOjir) — 2nqiqj j hi(r')hj(r + T , )^f t dT l . 


(IV.17) 
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Because hi(r') is zero for t' < 0, the lower limit of integration may be replaced by — oo. 
Since we are interested in the stationary cross correlation of a, and aj, the upper limit of 
integration may be replaced by oo, 


— OO — OO 


9 lwr • e lU}T gives 


Ra iaj (r) = 2Trq iqj J h^T^h^T + t')$ h dr' . 

— CO 

Substituting equation (IV.18) into equation (IV. 14) gives 

oo oo 

*«.•«» = J J hiir^hjir + T')e- iur dr'dT. 

Multiplying the integrand by t 

0*0 = J J hi(T l )e lwr ' hj(r + r^-^+^dr'dr. 

— oo — OO 

Integrating over r and then over t' gives 
where 

OO 

Hj { «)= J hjWe-^dt, 


(IV.18) 


(IV. 19) 


(IV. 20) 


(IV. 21) 


(IV. 22) 


is the complex frequency response function for the j th mode. H*{u) is the complex con- 
jugate of /fj(o;). Substituting equation (IV. 5) into equation (FV.22) and carrying out the 
integration give 

1 


*i(«) = 


u>j — ui 2 + 2iujj£jU> 


(IV. 23) 


Substituting equation (IV. 21) into equation (IV. 12) gives the power spectrum of the strain, 

(IV. 24) 


OO 00 


$ e€ (w) = $ /o £ ^ aiajqiqjH*(u)H j(u). 
i=l J=1 

In cases where is sufficiently small for all z, the double sum in equation (IV. 24) will be 
dominated by the terms in which i — j. In this case 




(IV. 25) 


1=1 


which is easily computed if the sum is truncated to a reasonable number of modes. 
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Section V . Simulation of the Random Time Series 

Given the power spectral density obtained analytically for a multi-degree- of-freedom 
system as described in the previous section, or obtained experimentally, one can construct 
a time series that will have the desired power spectrum. The method is primarily due 
to Rice [1] and it was modified to be computationally efficient by Yang as discussed by 
Shinozuka [2]. If $*(w) is the single sided power spectrum of the desired signal, s(<), then 
s(t) may be approximated by [1] 

N-i 

s(t) = y/2^2 [$i(w„)A u] 1/2 cos{u n t - <f> n ), (V.l) 

71 = 0 

where <f) n are uniformly distributed random numbers on the interval from 0 to 2n and 

u; n = nAuj, (V.2) 

with 

Au; = ui max /N } ( V.3) 

where u max is the maximum frequency in the power spectrum $*(a>), and N is the total 

number of terms in the summation in equation (V.l). 

Equation (V.l) simulates the time series as a distribution of sinusoidal signals having 
random phases. Unfortunately, this expression requires the computation of a large number 
of cosine functions at each desired value of the time, t. A considerable improvement in 
computational efficiency can be obtained by recasting equation (V.l) to allow the use of 
the Fast Fourier Transform. To accomplish this, note that equation (V.l) may be written 
as 

N - 1 

s(t) = Re {\/ 2 [^i(w„)Aw] 1/2 e ia, " t e - ^"], (V.4) 

n=0 

where Re[-] denotes the real part and i is 

If the simulated time series, s(t) is needed only at discrete values of time, t, then let 

si t = s(t k ) = s(kAt), (V.5) 

where the time duration between the equally spaced sample times is At. Evaluating 
equation (V.4) at t = t k gives 


N- 1 

s(kAt) = Re \\/ 2 ^ [$J(u>„)Aa;] 1/,2 e ,u '" iA^ e -, ^"], (V.6) 

n=0 


To satisfy the Nyquist sampling criterion, the time series, s(t) must be sampled at 
a high enough rate to obtain two samples during one period of the highest frequency 
component in the original input power spectrum, This gives 


At = 


27 r 1 

^m«i 2 


(V.7) 
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Substituting equations (V.3) and (V.7) into equation (V.6) gives 


N - 1 

s(kAt) = Re[V 2 Y, [$iK)A W ] 1 / 2 e-^"e 1 W i ], 

n=0 


(V.8) 


Equation (V.8) may be evaluated using a Fast Fourier Transform (FFT) algorithm by 
noting that given a discrete sequence, a n , the FFT provides an efficient means of computing 
Ak , where 


A k = Y ane- i27rkn ' N , 

71=0 


for k = 0,1,2,..., N- 1. 


(V.9) 


Equation (V.8) may be evaluated using a FFT by defining a sequence, 


°n — [^](w n )Au;] 1 / 2 e *^ n , for n < N — 1 
= 0, for n > N 


(V.10) 


Equation (V.8) may then be written as 

2JV — 1 

s k = Re[V2 Y a ne^], for k = 0, 1,2, . . . ,2N - 1, (V.ll) 

n=0 

where s k = s(kAt). Because we are taking the real part of the result of the summation, 
taking the complex conjugate of the right side of equation (V.ll) gives 

2N—1 

s k = Re[V2 Y for k = 0, 1,2, . . . ,2N - 1. (V.12) 

n — 0 

This is equivalent to 

s k = V2 Re[FFT(a n )], (V.13) 

where FFTf-] denotes the Fast Fourier Transform. Note that the length of the sequence 
a n is 2 N. ’ 
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Section VI. Comparison of Peak Counting and Rainflow Counting 

In this section time domain signals having a variety of spectral characteristics are 
processed to compare the fatigue life predicted using simple peak counting and rainflow 
counting. The goal is to identify features of the power spectrum that result in significant 
differences between the two damage counting methods. The methods described in the 
previous sections permit the construction of random time series corresponding to systems 
having any given spectral properties. From the discussion of section V one can construct 
the power spectrum of a resonant system having any number of resonant modes with 
specified resonant frequencies and modal damping ratios. It is also possible to create a 
random time series from a power spectrum that is obtained experimentally. 

Figures VI. 1 through VI. 6 show simulated time series and corresponding power spectra 
for systems modeled using the method of section V. These systems have nine resonant 
modes in the frequency range from 0 to 2000 Hertz. An attempt was made to create systems 
having spectral characteristics that mimic those seen in actual aerospace structures. The 
power spectra are shown in decibels by plotting 10Xo<7 10 (power spectrum). Both the 
actual input power spectrum and the power spectrum computed from the simulated time 
series are shown. Since the two curves are indistinguishable, it can be seen that the 
simulation procedure reproduces the desired spectral properties. Figure VI. 2 shows the 
power spectrum corresponding to the time signal in figure VI. 1, figure VI.4 contains the 
spectrum for the signal in figure VI.3, and figure VI.6 shows the spectrum for figure VI.5. 
As the power spectrum plots show, the time series differ only in the relative contribution 
of the first resonant mode. The figures showing the time series contain only the first 10 % 
of the simulated signal. 

The simulated time series were processed to estimate the fatigue life using peak count- 
ing and Rainflow counting. The ratio of the fatigue lives calculated by each method is 
shown on each plot. Between 4500 and 5000 damaging events were identified in each case 
shown. In these calculations the constant, c, in equation (III. 1.1) will cancel when the 
ratio of the fatigue lives are compared. The exponent, b in equation (III. 1.1) has been 
taken to be 6. 

In figures VI. 1 and VI.2 the first resonant mode has a relatively low amplitude. In this 
case the ratio of the fatigue lives computed using peak counting to Rainflow counting is 
.5910. In figures VI.3 and VI.4 the amplitude of the first mode is increased so that the time 
series shows a significant low frequency component. In this case the ratio of the estimated 
fatigue lives using the two methods is .1639. As the amplitude of the first mode is increased 
further as in figures VI.5 and VI.6, the ratio of the estimated fatigue lives reduces to .1333. 
The fatigue life predicted using Rainflow counting is then greater than that estimated by 
peak counting by roughly a factor of 7.5. When low frequency oscillations are present in 
the signal, peak counting significantly underestimates the fatigue life. 
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Figure VI. 1 Power Spectral Density 
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Figure VI.3 Power Spectral Density 
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Figure VL5 Power Spectral Density 
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VII* 1 Nonlinear Regression for Fatigue Parameter 

Any fatigue prediction procedure requires certain empirically derived constants to 
characterize the structure’s fatigue properties. In acoustic fatigue studies this data is 
usually obtained using a coupon test with narrowband random excitation. The S /N curve 
obtained from this test is then used to identify the fatigue constants. The application of this 
data to a broadband random load, as the structure is likely to experience in service, requires 
certain assumptions. The successful prediction requires that one is able to accurately 
account for the complex time histories encountered under broadband loading. The method 
of Rainflow counting as discussed in previous sections is intended to accomplish this. 

Although Rainflow counting has been used successfully in accounting for complex 
loadings it seems more direct (and therefore more accurate) to acquire the fatigue constants 
using a load history that more closely resembles the service environment. A procedure that 
could identify the necessary constants when complex loadings axe used could also be applied 
to determine fatigue properties from in-service failures. The goal of the present section is 
to develop a fatigue characterization method which can be applied with arbitrary loadings. 

The fatigue characterization scheme examined in the present effort begins with sets of 
data consisting of the time to failure along with the corresponding stress or strain power 
spectral densities. This data could be obtained either in tests of complex structures or in 
tests of simple coupon samples. Suppose that the measured fatigue lives are Tj and the 
corresponding spectra are $j(u>). Using the methods of reference [3] one can construct a 
simulated random time history and perform cycle counting to obtain a histogram, P ijt to 
describe the number of damaging events identified at each stress range, S, ; , for each power 
spectrum $y(cu). If the simulated time histories each have a duration t, then the damage 
experienced by the structure should be ° 


D i = t /Tj. (VII. 1) 

If there are M sets of fatigue lives and corresponding power spectra then the problem 
of estimating the fatigue coefficients, b and c, may be expressed as one of minimizing the 
total squared error, 

M / N C*. \ 2 

e(b, c) = Z(Dj-E ~T P 0 * AS ‘J (VII.2) 

j=l ' i=J C / 

where, b and c are constants depending on the materials, and b typically varies from 1 to 
6 while c from 10 20 to 10 30 . 

The problem of identifying b and c to minimize the error in equation (VII.2) may be 
simplified by replacing c by a new unknown, v, 

Co 

c = ^ (VII. 3) 

where, C 0 is a suitable constant. The choice of C 0 will be discussed later. 
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Substituting equation (VII.3) into equation (VII.2) gives, 


e(b, v) = J2 (dj - £ ■ A5, J 


or in general form, 


(**.**)- - E • A5 -) ; 


(VII.4) 


(VII.5) 


This is a modified mathematical model for identifying x\ and x 2 ( b and v) instead of b and 
c so that the property of the error function e(6, c) can be improved. 

To identify x x and x 2 , a combined Newton- Powell method is used in this project 
because of the special property of the error function e(x lt x 2 ) (or e(b,c )) which will be 
discussed in the coming section of this report. 

Given the nonlinear function e(xi,X 2 ), let 


X = (x 1 ,x 2 ) t = (VII.6) 

The Newton method solves this problem by expanding e(X) in a second order Taylor series 
about an estimate of the desired values of (x 1 ,x 2 ) T , X* = (x^x*) 71 . This procedure will 
enable us to obtain an improved estimate X fc+1 = (x* +1 ,x£ +1 ) T . In matrix form, 


e(X) = e(X*) + Ve T (X*) • (X - X*) + 1(X - X*) T ■ A(X*) • (X - X*) (VII.7) 

where, e(X*) is the functioned value of e at X*, and, 


V e (X fc ) 


de(xi 

dxi 

ae(x x ,x 2 ) 

dx2 


X=X k 


A(X‘) 


a z e(n r 2 ) 


dxi 


9 2 e(xi,r 2 ) 

dzidxi 


8 2 e(xi,x 2 ) \ 
dxidx^ I 
a 2 e(z 1 ,Xj) I 
/ 


x=x k 


and, 


5e(x!,x 2 ) 

dxi 




j = i 


t=i 




/n(x 2 • Sjj) • (x 2 • 5,j) ri 
Co 


Pii 



(VII.8) 


(VII.9) 


(VII.ll) 
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(VII. 12) 


de(x u x 2 ) 

dx? 


M 


N 


= 2 £(^-£ {fi #^o'As) 


- « c ° 

" ( J 1 ' Sjj) • (x 2 ' Siy)* 1 " 1 

Co 


-E 


d 2 e(x 1 ,x 2 ) 


dx\ 


«= 1 

Mr, N 

= 2£ (*>-£ 


P tJ • AS,) 


j=i L 


i=l 


(f2 • Sr,') 11 

Co 


5 0 • AS;) 


( _ Y' [Mf2 • S^)] 2 • (a^2 • S,j) Tl 
Co 

Mf 2 • S',;) • (x 2 • S.j)* 1 


E 

i=i 
N 

+ (-£ 

1=1 


3 ij • AS,) , 
P.i • AS,) 


^ (*! • Sjj) • (x 2 • $,■)»»-* 


dx\dx 2 


1 V 1 r y iV 

= ^£ (-£ 

jsl X 1=1 


Co 


TV 


ln(x 2 ■ 5jj) - (;r 2 - 5V;) Zl 


£ 

i=l 


• AS.) 
P„ ■ AS;) 


JV 


+ *>-£ 


1=1 


(x 2 • S,j) ri 
Co 


P.j ■ A S;") 


£ 

1 = 1 


[Sjj + X! • /n( x 2 • Sij)] ■ (x 2 • SijY 1 1 
C 0 


Pij • AS, 


d 2 e(xx,x 2 ) 
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To find the minimum point, let 

V e (X) = A(X k ) • (X - X*) + ye(X*) = 0, 


(VII. 13) 


, (VII. 14) 


(VII. 15) 


(VII.16) 
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Solving equation (VII. 16) for X leads to an iterative procedure to determine the next 
estimate X* +1 , 

X = X* +1 = X* - A(X*) -1 • V e(X fc ), (VII. 17) 

where, A(X*) -1 indicates the inverse of matrix A(X*). It is assumed that A(X*) is a 
positive definite matrix and its inverse exists. 

As it will be discussed in the coming section, for the error function given in equation 
(VII.5), the matrix A(X*) will become singular when X* is near the minimum point of 
e(X). This suggests that, although the Newton method does an excellent job when X* 
is far away from the minimum point for this problem, it will fail when X* is near the 
optimum approximation of b and c. To obtain the optimum estimate of 6 and c, we use 
the Powell method [4] to find the optimum values of b and c when the Newton method 
begins to fail to do so. 

Here is how Powell method goes: 

Initialize the set of directions S' to the base vectors, 


S 1 = (1,0,0, ...,0) T , s 2 = (0, 1 , 0 , ..., 0 )^ 


s N = (0, 0 , 0 , ..., 1) J 


Now repeat the following sequence of steps (“basic procedure”) until the function stops 
decreasing: 

• Save the initial estimation of X as X°. 

• For i = 1 , TV, move X’ -1 to the minimum along direction S' and call this point 
X*. 

• For i = l,...,jV- 1, set S* «- S i+1 . 

• Set S N <- X N -X°. 

• Move X N to the minimum along direction and call this point X°. 

In addition, some modifications have been added to the basic procedure to ensure a 
reasonable rate of convergence. 


VII.2 Identifiability of Parameters b and c 

As mentioned in the previous section, the error function e(6, c) possesses some special 
properties. This affects the identifiability of parameters b and c. First, we can make a 
simple observation. When carrying out a fatigue test, we excite a structure by using a 
random loading with a specified stress power spectrum, such as what was done in Cleven- 
son’s [5] and Ramesh’s [6] fatigue tests in which each specimen was loaded at a constant 
root-mean- square stress level with a given shape of stress spectrum (as shown in figures 
(VII. 1) through (VII.9)) until rupture occurred. In order to estimate the fatigue parame- 
ters b and c using the method developed in this project, the data of the root-mean-square 
stress and the stress spectral shape are used to construct the stress spectrum under this 
certain test condition. We will assume that a set of fatigue lives, Tj,j = 1,2 have 
been measured and that the corresponding stress power spectra, Sj(v),j = 1,2 
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axe known. The methods of reference [3] are used to construct a simulated random time 
history during a record of duration r (as shown in figures (VII.10) through (VII.18)) for 
each fatigue test. Rainflow cycle counting is then used to obtain a histogram of the levels 
of damaging cycles, P\j (as shown in figures (VII. 19) through (VII. 27)), where the index i 
corresponds to the damage level and j indicates the test specimen number. The value of 
Pij ■ A Si, as in equation (VII. 2)), describes the number of damaging events identified at 
each stress range during a record of duration r. Therefore, b and c should be constants 
such that, 

T N S b 

D > = t’ = T, -T p 0 • j = (VII.18) 

i = 1 c 

Since the fatigue characteristics of a given material are assumed to be described by the 
two constants, b and c, it is hoped that for any given set of test data, a unique pair of values 
of b and c may be determined. Unfortunately, due to the nature of fatigue measurements, 
an error in the determination of b can be partially off-set by adjusting c. This problem 
comes about when b and c are to be identified in tests conducted with simple deterministic 
load histories as well as in the more complicated situation considered here. 

First consider the case where the load history consists of a pure tone so that the number 
of damaging cycles is easily determined. A typical S/N curve for sinusoidal loading [5] is 
shown in figure (VII.28). 

Figure (VII.29) shows estimated fatigue curves having significantly different values of 
b and c. The values of b and c for each curve are shown in the figure. It is found that if one 
parameter, either b or c, is chosen incorrectly, then the other parameter may be adjusted 
to compensate for much of the error. Figure (VII.29) shows that values of c varying over 
two orders of magnitude can give reasonable results as long as 6 is adjusted accordingly. 

A similar affect to that seen in figure (VII.29) is also observed in the more general 
loading case considered here. To illustrate the c ne ar 1 dependence of b and c in this case, 
consider the problem of finding b and c when there is only one set of measured fatigue test 
data. In this case equation (VII.2) becomes, 

e(6, c) = ~ 'fPij 1 AS,) , j = 1- (VII.19) 

Figure (VII.30) shows the error, e(6,c), as computed using equation( VII.19) for a 
range of values of b and c. The figure shows that for each value of 6, a value of c may 
be found to cause the error to be zero. One can then use equation (VII.19) to construct 
a relationship between b and c which gives e(b,c) = 0, as shown in figure (VII.31). This 
relation between b and c may be obtained for each fatigue test specimen if a series of tests 
are conducted. Figure (VII. 32) shows the relationships between b and c which give zero 
error for 4 different tests [5] having spectra identified as A, B, C and D. It would be ideal 
if all of the lines intersected at a particular combination of b and c. The figure shows that 
the lines indicating a relationship between b and c corresponding to the different fatigue 
tests are almost parallel. There is thus no ideal combination of 6 and c for this data. This 
results in the error function e(6, c) being illconditioned. 
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Another way to illustrate the near dependence of b and c in typical test data is to 
plot the error as a function of b and c as in figure (VII.33). The minimum of e(b, c) is seen 
to consist of a curve rather than a point. It is obvious that there exists a narrow region 
in which the value of error function is almost kept unchanged. The minimum values of 
the error function with respect to b and c are shown in figure (VII.34) and figure (VII.35), 
respectively, in which it can be found that the minimum value of e(6, c) just changes very 
slightly while b varies from 3.5 to 5.5 and c from 10 19 to 10 28 . This results in the failure of 
Newton Method for this problem when the iteration point approaches this narrow region 
because where the matrix A(X*) will become near singular. This also suggests that a 
range of values of b and c will give almost equally good estimates to the data. 

There are many situation where we can find some general empirical relationship be- 
tween fatigue constants [7], When the S/N curves for some steel alloys are plotted in 
nondimensional form using the endurance limit, S e , and the ultimate strength, S u , they 
tend to follow the same curve shown in figure (VII.36). A power relationship is recom- 
mended to be used to estimate the S/N curve for steel: 

S = 10 c • A -6 (for 10 3 < N < 10 6 ) (VII.20) 

where the exponents, C and b, of the S/N curve are determined using the two defined 
point shown in figure (VII. 3 7): 

L — 1 7 ■S’iooo ~ , (•S'looo) 2 

o=—l°9io—^—, C = log 10 — - (VII.21) 

where, alternating stress level corresponding to a life of 1000 cycles, S 1000 , can be de- 
termined as 0.9 times the ultimate strength, and when the estimate for S e is made as 
S e = 0.5S U , the S/N curve is defined as 

S = 1.625 U JV“ 0 085 (VII.22) 

This suggests that for certain kinds of steel, the S/N curve can be approximately defined 
as only depending on the ultimate strength S u , which also suggests that fatigue constants 
are not strongly independent. This results in the ill-conditioning of the error function. 

In order to overcome the ill-conditioning, it is necessary to make a dimensional trans- 
formation of variables. When using equation (VII.2) as the mathematical model for this 
problem, the contours of the error function (with b and c as two parameters) near the min- 
imum point are often curved as in figure (VII. 38) ( b and c are both in linear scale because 
they behave like th is when we estim ate the parameters using Newton Method or Powell 
Method) . Here the projection of the contour has substantial length. After a nonlinear 
dimensional transformation 

c ^ (vn.3) 

is made as in Section VII. 1, the contours of e(b,v) near the minimum point of the error 
function are a s shown in figure (VII. 39). This makes the estimation of b and c much more 
stable than that before the transformation. It has been found that a suitable value of the 
constant, Co, is Cq = 2xl0 22 in this case. 
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VII.3 Application to Existing Fatigue Data and Discussion 


In order to show the advantage of the method developed in this project, it has been 
applied to some existing measured fatigue data published in both references [5] and [6], in 
which some statistical properties of the random loading are used to estimate the fatigue 
constants and to predict the fatigue lives of the structures studied. 

Fatigue tests of aluminum- alloy specimens under stationary and Gaussian random 
loadings having a zero mean value were conducted in reference [5]. The random loadings 
were described by power spectral densities of the form $(u>) = -Kiu> -2 , J£ 2 u>°, K 3 u> 2 in 
a certain frequency range, where K is an arbitrary constant ( typical stress spectra are 
spectrum A, B, C, and D as shown in figure (Vll.l) through (VII.4)) and the fatigue 
lives were determined for statistical parameters of the random loadings such as root-mean- 
square nominal applied stress, power spectral shape, mean number of zero crossings per 
unit time, and mean number of peak loads per unit time, etc. The corresponding fatigue 
data are listed in table (VII. 1), which shows root-mean-square nominal applied stress, 
Srmai * n P s h fatigue lives , tf, in seconds, a, the ratio of zero crossing rate to peak rate as 
defined later in equation (VII. 27), and the corresponding stress spectrum used in the test. 

The measured data in reference [6] were obtained by fatigue tests of a number of 
structural steel specimens. Five typical power spectral shapes were selected to generate 
various random loading wave forms and the idealized shapes are spectrum A, B, C, D, 
and E as shown in figures (VII.5) through (VII.9) by setting the filters in different ways. 
Fatigue tests were conducted under different power spectral shapes, measured on the basis 
of central frequencies, positive zero crossing frequencies, and frequencies of maxima. Table 
(VII.2) shows the corresponding fatigue data which gives power spectral shapes used, root- 
mean-square applied stress, S rm3 , the fatigue lives, tf, and the value of a. 

To apply the method developed in this study, time histories for each power spec- 
trum mentioned above are generated (as shown in figures (VII.10) through (VII.18)) and 
a Rainflow counting scheme is used to calculate the damage accumulation in each case 
investigated. Some typical plots of damaging densities of the stress ranges identified by 
Rainflow counting are shown in figures (VII. 19) through (VII.27). 

The simulated damage probability densities and corresponding measured fatigue lives 
may be processed as described in equations (VII. 1) through (VII.5). By applying the 
scheme discussed in the previous section to some of the measured fatigue data in both 
references [5] and [6], fatigue constants b and c Eire identified and then used to estimate 
the fatigue lives for all other test specimens. Figures (VII.40) and (VII.41) indicate the 
predicted fatigue lives using the identified fatigue constants, which show that, with this 
method, the identification of fatigue constants and the prediction for fatigue lives based 
upon the Rainflow counting technique can be made with reasonable precision. 

Some comparison are also made to show the efficiency of this method: 

First, the effect of power spectral band width is considered. Let a be the ratio of zero 
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crossing rate to peak rate. This may be computed if the power spectrum, $(u>), is known, 


a = 


r« 2 $(u)(L 


(j™ w$(u)duj uj 4 $(u)du?J 


(VII.27) 


The a values for tests conducted in both references [5] and [6] vary from 0.586 to 
0.992, as shown in tables (VII. 1) and (VII.2). This indicates that the data corresponds 
to broadband and narrowband random loadings. A power spectrum with a=0.586 is con- 
sidered to be a broadband random loading and when or = 1, the loading is narrowband. 
The difference between narrowband and broadband spectra can be seen from typical time 
histories shown in figures (VII.12) and (VII.14), which have a values of 0.695 and 0.992, re- 
spectively, and the corresponding spectra are shown in figure (VII.42), where both spectra 
have the same unit root-mean-square stress value. 

The application of the regression method to each power spectrum used in the tests is 
investigated. From table (VII. 1), 32 sets of fatigue data are shown which are taken from 
reference [5]. The tests were conducted using four different spectrum shapes identified as 
A, B, C, and D. For spectrum shapes A, B, C, and D there are twelve, four, ten, and six 
measured fatigue lives, respectively. 

Figures (VII.43) through (VII.46) show the different cases of fatigue life prediction. In 
each of these figures, data from one of the four spectrum shapes was used to identify b and c. 
The data used to perform the regression are shown as circles. In each figure, predictions of 
the fatigue lives for the data of the remaining three spectrum shapes are shown as squares. 
Figure (VII.43) shows twelve circles, representing a comparison of twelve predicted and 
measured fatigue lives obtained under spectrum A. These twelve circles indicate the fatigue 
test data used to identify b and c while the squares indicate predicted fatigue lives based 
on these constants for the rest of the fatigue data listed in table (VII.l). Figures (VII.44), 
(VII. 45) and (VII.46) show predicted lives in the same way as in figure (VII.43), in which 
fatigue constants are identified from fatigue data obtained under the random loadings with 
spectra B, C, and D as shown in figures (VII.2), (VII. 3), and (VII. 4), respectively. 

We can also find the similarity in predicting fatigue lives using identified fatigue con- 
stants from figures (VII.47) through (VII. 51) in which, for fatigue tests conducted in refer- 
ence [6], fatigue constants b and c are identified from those fatigue test data under random 
loadings with five different stress spectrum shapes A, B, C, D, and E, respectively. It can 
be found from figures (VII.43) through (VII.51) that, with this method, the prediction 
for fatigue lives under the random loading with certain power spectrum can be made al- 
most equally well by using fatigue constants b and c estimated from fatigue tests under the 
loading with very different power spectral shapes even under narrow band random loading. 
Figures (VII.43) and (VII.47) show the extreme cases in which fatigue lives are predicted 
using fatigue constants identified from the test data measured under narrow band random 
loadings with at values of 0.92 and 0.99, respectively. It means if the fatigue constants b 
and c are determined using narrowband test data, the prediction can also work well for 
fatigue lives under broadband random loading. 
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Secondly, since Rainflow counting considers the effect of each stress level, it seems 
to be possible to identify fatigue constants b and c by using only one or a few more sets 
of fatigue test data instead of large amount of fatigue tests. Figures (VII.52) through 
(VII. 60) show the predicted fatigue lives by using only one set of fatigue test data chosen 
for each type of spectrum. The results show very good agreement. It is expected that, on 
the consideration of economy and for the easy implementation of the fatigue test, one can 
conduct a fatigue test using only one spectrum shape which is specially designed so that 
damaging events occurs more evenly at each stress level. 

The practical advantage of using Rainflow cycle counting to calculate the damage 
accumulation in predicting fatigue life can be pointed out here. By using the Rainflow 
counting technique, one can consider the effect of each stress level on fatigue life in both 
identification of fatigue constants and prediction of fatigue lives even if only one set of test 
data is available. Figure (VII. 5) indicates a narrow band spectrum shape corresponding 
to a time histoiy as shown in figure (VII. 14) which can be found to be a typical narrow 
band random time history. But, from figure (VII.23), we can find that it includes as much 
information about damage events at each stress level as that in the broad band case. This 
is one reason that we can get almost equally good predictions of fatigue life by using fatigue 
constants estimated from fatigue tests under random loadings with very different power 
spectral shapes. 


VII.4 Conclusions 

The common method of applying narrowband random S/N data is to simply consider 
the structure under study to also respond in a narrowband, random fashion. In practice, 
however, a structure is expected to encounter a more complicated broadband random 
excitation. A primary goal of this effort has been to develop a method of identifying the 
fatigue constants when a realistic broadband random load is applied. It was felt that this 
would lead to more accurate predictions. We have found that when Rainflow counting 
is used to identify damaging events and when the fatigue prediction is performed using 
the methods outlined above, one can make reasonably accurate predictions of broadband 
fatigue life using narrowband fatigue data. 

Because the success of the present method depends strongly on the accurate iden- 
tification of damaging events in a complex time history, one can conclude that Rainflow 
counting can be used reliably. The regression procedure we have described can be useful in 
situations where narrowband fatigue data is not available such as when fatigue constants 
axe to be determined from in-service data. 
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Figure Vn.7 Power Spectrum C [6] 




Figure VII. 10 Time History of Spectrum A [5] 










Figure VII. 20 Damaging Probalility Density 
of Spectrum B [5] 



Figure VII. 21 Damaging Probalility Density 
of Spectrum C [5] 



Figure VII. 22 Damaging Probalility Density 
of Spectrum D [5] 
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Figure VII. 24 Damaging Probalility Density 
of Spectrum B [6] 



Figure VII.25 Damaging Probalility Density 
of Spectrum C [6] 



Figure VII.26 Damaging Probalility Density 
of Spectrum D [6] 
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Figure VII. 27 Damaging Probalility Density 
of Spectrum E [6] 
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Sinusoidal Loading [5] 
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Figure VII.30 Error Functional Values vs. b and c 



Figure VII. 31 Relationship between b and c 



Figure VII.32 Zero Error for 4 Different Tests [5] 



Figure VII.33 Error Functional Values vs. b and 
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Figure VII. 37 Generalized S/N Curve for Wrought 
Steels on Log-Log Plot [7] 
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VTI.40 Measured and Predicted Fatigue Lives [5] 



Figure VII. 41 Measured and Predicted Fatigue Lives [6] 
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Figure Vn.42 Samples of Broadband and Narrowband Spectra 



Figure VII. 43 Measured and Predicted Fatigue Lives [5] 
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VII. 46 Measured and Predicted Fatigue Lives [5] 



Figure VII.47 Measured and Predicted Fatigue Lives [6] 




Figure VII.48 Measured and Predicted Fatigue Lives [6] 



Figure VII.49 Measured and Predicted Fatigue Lives [6] 
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Figure YII.51 Measured and Predicted Fatigue Lives [6] 
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VII. 52 Measured and Predicted Fatigue Lives [5] 



Figure VII. 53 Measured and Predicted Fatigue Lives [5] 
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Figure VII.54 Measured and Predicted Fatigue Lives [5] 
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Figure VII. 55 Measured and Predicted Fatigue Lives [5] 
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VII.56 Measured and Predicted Fatigue Lives [6] 



Figure VII.57 Measured and Predicted Fatigue Lives [6] 



Figure VII.58 Measured and Predicted Fatigue Lives [6] 



Figure VII.59 Measured and Predicted Fatigue Lives [6] 
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Figure VII.60 Measured and Predicted Fatigue Lives [6] 


Table VII.l Fatigue Data for Random Loading [5] 















































Table VH.l - Continued [5] 

































































Table VII. 2 Fatigue Data for Random Loading [6] 


SPECTRUM APPLIED STRESS FATIGUE LIVES ZERO CROSSING RATE/ 
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13924 



14940 
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16680 

A 

17840 

B 

14940 

B 

14940 

C 

14940 

D 

14940 

D 

16100 

D 

18420 


tp sec. 

PEAK RATE 

37394 

0.992 


20170 0.992 

10522 0.991 


12924 

25680 

26330 

24320 

53870 

29390 

13770 
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Section VIII 


Powell’s Method 


Powell’s method has the characteristic that when it is applied to a quadratic form, 
it choose conjugate directions (defined in the following paragraph) in which to move. 
Hence, the rate of convergence is fast when the method is used to minimize a general 
function. In addition, Powell has added a modification to the basic procedure to ensure a 
reasonable rate of convergence when the initial approximation is quite poor. 

Conjugate Directions 

Consider the problem of finding the minimum value of the quadratic function 

/(X) = Jx t AX + B t X + C, (VIII.l) 

where, A is positive definite and symmetric. In two dimensions, the level curves /(X) = K, 
for different values of K, are concentric ellipses as shown in figure (VIII.l). Suppose that 
we search for a minimum from the point X° in the direction S 1 , that this minimum occurs 
at X 1 and that C is the optimal point. We say that the direction X 1 — C is conjugate 
to the direction S 1 since, for any ellipse /(X) = K, the diameter through X 1 is conjugate 
(in the usual geometrical sense) to the diameter parallel to S 1 . The idea of conjugate 
directions is easily extended to n dimensions. 

Let a set of non-zero search direction vectors S^S 2 ,..., S n be conjugate to a given 
positive definite matrix A, if the following condition is satisfied, 

S ,T • A • S' = 0, for all i^j, (VIII.2) 

where the matrix A is as defined in equation (VIII.l). 

We can reach the minimum in at most n exact one-dimensional searches, if the i-th es- 
timate, X‘, is generated by the initial guess, X°, and the conjugate directions S 1 , S 2 , ..., S'. 
The estimates of X at each iteration may be written in the form, 


X 2 = X 1 + a 2 S 2 = X° + 01 S 1 + o 2 S 2 , (VIII. 3) 

X 1 = X° +a x S l , 

i 

X 1 ' = X*" 1 + a; S’ = X° -(- 5^ a > si , 

j= i 

where oti are optimum step size it should go from X' -1 to X‘ along direction S*. 

A search technique which finds the minimum of /(X) defined by equation (VIII.l) 
in a finite number of one-dimensional searches will have considerable merit for a general 
optimization problem, since it should converge rapidly once the optimum is approached. 
The basis for these methods can be seen in the geometric construction shown in figure 
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(VIII. 1). Here any point, X°, is taken as the initial base point, and a search is made in an 
arbitrary direction (shown here as S 1 ) until a minimum, X 1 , is found. A second point, X 2 , 
not on the line through X° and X 1 , is now taken and a search carried out for a minimum, 
X 3 , on the line parallel to that through X° and X 1 (S 1 ). If the contours are ellipsoidal' 
as defined by equation (VIII. 1), a line connecting points X 1 and X 3 constructing a new 
direction, S will pass exactly through the true minimum. A one- dimensional search 
along this connecting line will then reveal the desired optimum. We can therefore find 
the optimum without further iteration, provided the objective function has the quadratic 
form given by equation (VIII. 1). If the objective function is not quadratic, then of course 
convergence is not assured in a finite number of iterations, or even after a very large number 
of tries. Since most problems of interest are not expressed in terms of quadratic objective 
functions, we might inquire as to the value of quadratic convergence. Most functions, 
however, approach quadratic behavior in the neighborhood of their optimum, and since 
quadratic functions can have very elongated elliptical contours, quadratically convergent 
methods are capable of coping with such shapes. In essence, it seems very reasonable 
to expect a quadratically convergent procedure to be effective in searching nonquadratic 
object functions. 

Powell's Basic Method 

Initial directions in Powell’s method are usually taken to be the coordinate directions, 
S 1 =(1,0,0,. ,.,0) T , S 2 = (0,1,0, ...,0) T , ..., S n = (0,0,0, ...,l) r (VIII.4) 

then repeat the following sequence of step until the function stops decreasing: 

(1) Save the initial estimation of X as X°. 

(2) For i = 1, ..., n, move X ,_1 to the minimum along direction S‘ and call this point 
X‘. 

(3) For i = 1, ..., n - 1, set S‘ S ,+1 . 

(4) Set S n ^ X n - X°. 

(5) Move X n to the minimum along direction S n and call this point X°. 

But, if one uses Powell’s method to solve 

/(X) = x 1 2 -z 1 x 2 + 3x 2 2 (VIII.5) 

starting at the point X° = [1,2] with S 1 = [1,0], and S 2 = [0, 1], then it will be found that 
when searching from X° in the direction S 1 for a minimum, it turns to be the point X° 
itself. Thus X = X® and the method breaks down: the directions of search in subsequent 
iterations are restricted to a subspace which does not contain the direction S 1 . The given 
problem can not therefore be solved by this method. 

Powell's Method 

Powell modified his basic method to overcome the type of difficulty encountered above 
by allowing a direction other than S 1 to be discarded after each iteration. In this way, the 
n directions of search can be chosen so as to be always linearly independent; in some cases, 
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the same n directions are used for two successive iterations. For each iteration starting 
with initial point X and n directions S-^, j = 1,2, ...,n, it proceeds as follows: 


(1) For i — 1,2,..., n, search from X* 1 in the direction S' for a minimum at X*. 

(2) Find 


A = - f( x ’) 

(VIII.6) 

i.e. q is the value of i which has a maximum value of A. 


(3) Define 


fi = /(X°) h = /(X") 

(VIII.7) 

and evaluate 


h = /( 2X n - X°) 

(VIII.8) 

(4) If either 


V '*S 

w 

IV 

> 

(VIII.9) 

and/or 


2(/i - 2 h + h)U i -h- A) 2 > A (h - f z )\ 

(VIII. 10) 

use the old directions S', i = 1,2, ...,n for the next iteration and use X n 
point of the next iteration, X°. Otherwise, use rule (5). 

as the initial 

(5) Set 8 - X" - X°; move X n to the minimum along the direction 8 and call this 
point X°; and take the following sequence as the directions for the next iteration 

S 1 , ..., S 7-1 , S g+1 , ..., S", 8. 

(VIII.11) 


In (4), the first inequality considers the value of a tentative value of the objective 
function (error function) fa found by exploring ahead along the new search direction. 
Here Powell tests the objective function at the point X = 2X n - X°, which is a point 
along the direction S n , located the same distance from X" as the original base point X°. 
If this tentative point in the new direction does not produce a better value of the objective 
function than the original base point, the new direction is rejected without further tests 
as being a poor prospect. 

The second inequality is used to determine if the function may not be rising sharply 
( when seeking a minimum ) in the move from the point X n to the tentative point X 
after dropping sharply from X° to X n . In other words, the test is used to determine if 
the direction S n is pointing across a deep valley, and if it is, we reject it as a reasonable 
direction of search, retaining the previous search directions for use in the new stage. 


Reference 

[1] L. Cooper and D. Steinberg 1970, Introduction to Methods of Optimization W. 
B. Saunders Co. 
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IX. 1 Governing Equation for a Nonlinear Beam - Rainflow Counting 


The results presented in the following are part of an investigation of the effects of 
nonlinear vibration on the fatigue of structures. Nonlinear effects are important in fatigue 
studies because it is usually large deflection response that causes damage and nonlinearities 
are more prevalent when the deflections are large. While most structural systems are not 
subjected to sufficient excitations to elicit nonlinear response in service, it is not uncommon 
to perform experimental studies of the fatigue life with artificially high vibration levels. 
By testing at high vibration levels the rate of damage accumulation may be accelerated so 
that the test can simulate a long service life (say 10 to 20 years) during a normal 8 hour 
working day. In order to properly use the results of the accelerated test, however, it is 
important to know what influence nonlinearities may have on the results. The goals of the 
present investigation are to examine the impact of common nonlinear effects on fatigue life 
and to develop practical tools for including nonlinearities in fatigue predictions. 

In the present section the method of Rainflow' counting is applied to estimate the 
fatigue life of a nonlinear beam. The model of the response of the beam comprises an 
extension of the linear beam analysis presented in section II.2 to account for nonlinearities. 
The equation of motion for the deflection, W , of a base excited, elastic beam including the 
effect of nonlinear stretching strain at large deflection is 


pAW(x,t) + 
Ebh 


Ebh 3 . 


■W x 


l 


Wr 


12 

f W 2 x {t,t)dt + r ] W = -pAW 0 . 
Jo 


(IX.1.1) 


where p is the density, A is the constant cross sectional area, E is Young’s modulus, b 
is the width, h is the thickness, / is the length, and r/ is a viscous damping coefficient. 
The subscript x denotes partial differentiation with respect to x and dots denote partial 
differentiation with respect to time, t. W is measured relative to the displacement of the 
beam ends which are assumed to have equal prescribed acceleration \V 0 , and have zero 
rotation. W a will be assumed to be stationary Gaussian white noise. In the derivation of 
equation (IX.1.1) it is assumed that there is no axial displacement at the beam ends. The 
derivation is presented in section IX.4. 

The relation between the strain and displacement of the beam may be shown to be 


1 




dx 2 ’ 


(IX. 1.2) 


where z is the vertical distance from the mid-line of the beam. Since the boundaries 
prevent rotation at the beam ends, the maximum strain will occur at x = 0, and x = / 
and at the extreme fibers at z = ±h/ 2. The maximum strain is then 


where x = 0 or x = /. 



hd 2 w(x,t) 

2 dx 2 ' 


(IX.1.3) 
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An approximate solution of equation (IX. 1.1) may be developed by expanding the 
displacement response in a finite set of the eigenfunctions for a linear cl am ped- clamped 
beam, 

N 

W(x,t) = ^ai(t)4>i(x), (IX. 1.4) 

«=1 

where N is a finite integer and a,(t) are unknown functions of time. For the linear beam, 
the eigenfunctions are given by 


4>i(x) = cos (pix/1) - cosh (p,x//)+ 
£>,(sin(pjx//) — sinh (p,x//)) 

(IX. 1.5) 

where p, and D t are given in Table II. 2.1 for the first six eigenfunctions. 


Substitution of equation (IX. 1.4) into (IX. 1.1) leads to a coupled system of nonlinear 
ordinary differential equations for the a,(<), 

N 

a, + Cdi + ujfai + Q(a ) Iij a j = F \W 0 {t ), * = 1, 2, . . . N, 

i=i 

(IX. 1.6) 

where a is a vector containing the a ,•(<), and 


F N N 

= XX 

P «=1 jm 1 

(IX. 1.7) 

II 

| _ 

(IX.1.8) 

2= £fc 2 /PA< 

' 12 p \ 1 ) ’ 

(IX. 1.9) 

T _ 1 [‘ dM*) dtjfx) , 
l] 0 dx dx ix ' 

and 

1 t l 

Fi = -- J <j>i(x)dx. 

(IX.1.10) 

(IX.1.11) 

Note that due to the symmetry of the beam F t is zero if i is even. As 
steady-state response corresponding to even numbered modes is zero, 

a result, the 

lim Off (f ) = 0 if £ is even. 

t — >00 

(IX. 1.12) 

Once the arj(t) are determined, the displacement, W (x, t), is obtained from equa- 
tion (IX. 1.4) and the maximum strain may be calculated using equation (IX.1.3). 


107 


IX. 2 Quasi-Simulation for Multi-Mode Nonlinear Systems 

In this section, we shall present the method of equivalent linearization for solving 
equation (IX. 1.6). A simulation procedure of the resulting equivalent linear system is then 
developed to generate the maximum strain time history for the purpose of fatigue analysis. 

The method of equivalent linearization consists of approximating equation (IX. 1.6) by 
an equivalent linear system, 


N 

<*i + C<*. + ^2 Ki i ai + e >te) = F iWo(t), 1 = 1,2, ... N, (IX. 2.1) 

i - 1 

where the constants K t] are chosen to minimize the error, e,(o), in some sense. PYom 
equations (IX. 1.6) and (IX.2.1), the error term e,(a) can be written in a vector form 


N N 

e = {e,(o0} = <j J2 K ii a i ~ ~ <?(«) £ [ ■ 

>=i j=i 


(IX.2.2) 


Minimization of the steady-state mean square value of the e with respect to Kij leads to 

5E [e T e] 


dKn 


= 0. 


(IX. 2. 3) 


Expanding equation (IX.2.3), and using equation (IX.1.7) give 


N 

^ ^ = E[orjCKj] + 

m=l 

E N N N 

2 ~ ^2 ^2 ^2 

F m= 1 k=l 1=1 


(IX.2.4) 


The solution for the o : ’s from equation (IX.2.1) with the error term neglected is known to 
be jointly Gaussian. The fourth order moments in equation (IX.2.4) can thus be expressed 
in terms of lower order moments, 


EtajOfjnafcaj] = E[a ; o: m ]E[afcar/]-|- 

E[<*><**]E[a m a!/] + E[ajO;j]E[a m afc]. 


In deriving equation (IX. 2. 5), we have used the fact that all the odd order steady-state 
moments of the a,’s are zero. Substituting equation (IX.2.5) in equation (IX.2.4) gives 

N 

^ 1 -^imE[a m aj] = w^E[ajaj] + — x 
m=l 2 P 

N N N (IX. 2. 6) 

EXE Iim A/{E[ofja m ]E[a'jta/] + 2E[aja*]E[a m a/]}. 

T7l=l k=l 1= 1 
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Equation (IX. 2. 6) can be used to determine the components of Kij. However, one has 
to solve simultaneously for E[q,Qj] and Kij using an iterative procedure as described in 
the following. To simplify the discussion, we introduce the following matrix notations: 


R = {E[a,-orj]} , K = {Kij} , I ={/,,}, 


n 2 = 


/“> f 


V 


U!n 


\ 


"I' 


Equation (IX. 2. 6) may then be written in a compact form 


(IX. 2. 7) 


K • R = f} 2 • R + — {(I: R)I • R -f 2(1 • R) • (I • R)} , 

Ip 


(IX. 2. 8) 


where 


N 


N N 


( K ■ R)o = X *<»*»>. (* R ) = X E 

m=l «=1 j=i 

Solving equation (DC.2.8) for K, we have 


>2 E 


K = SI 2 + -{(I:R)I + 2(I.R) ■!} 


(IX. 2. 9) 


Notice that in general, K is not diagonal. Neglecting the error term, e;(o), in equation 
(IX.2.1), then leads to a coupled system of linear ordinary differential equations. With 
e i(o) neglected, equation (IX.2.1) may be solved using modal analysis. 

Consider the equation 

a + K« = 0. (IX.2.10) 

Let yf (* = 1, 2, . . . , N) be the eigenvalues of the matrix K in equation (IX.2.10) and let 
$ be the normalized matrix of eigenvectors such that 


(l\ 


$ T • $ = J, $ r K -3> = 


7f 


where X is the unit matrix. Also let 


\ 



(IX. 2. 11) 


a = $/3, f=$ r F. (IX. 2. 12) 

If the error term, e,-(a), is neglected, equations (IX.2.1) and (IX.2.12) lead to the equations 
for modal coordinates /?, 


Pi + CPi + 7 , 2 Pi = fiW 0 , i = 1, 2, . . . N. 


(IX.2.13) 
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This uncoupled system of equations will be used in the following to simulate the time 
domain response of the nonlinear beam. Define a matrix of second order moments of the 
modal coordinates, 

Tl = {E [fiifo]} , (DC. 2. 14) 

where the steady-state moment E [0i0j] can be obtained as 




(7 l ? -7j) 2 + 2C 2 (7, 2 + 7?)’ 


(IX.2.15) 


and where is the single sided power spectral density of W 0 with units of acceleration 
squared per Hertz. The matrices R and 1 Z are related by 

R = $ K & T . (IX. 2. 16) 


Using equations (IX. 2.9) through (IX.2.16), the iterative procedure for determining 
R and K may be stated as follows: 

Step 1. (Zero order approximation) Assume that 

EfajOj] = 0, K, j = 0 for i j. 

Then, from equation (IX. 2. 6), 


K “ = "? + + 2^)E[ a ;i 

H i- i 


(IX.2.17) 


By neglecting the error term, 
elements of K tJ , we obtain 


e '(^-) * n equation (IX. 2.1) and neglecting the off-diagonal 


Ku = 


4CE[af] ' 


(IX.2.18) 


Equations (IX.2.17) and (IX.2.18) can be readily solved iteratively by starting with an 
initial value of ha = w?. Equation (18) may then be used to determine E[nJ] and the 
result substituted into equation (IX.2.17). Evaluation of (IX.2.17) then produces an up 
dated value of Ku. This process may be repeated until K tt and E[o?] converge. 

Step 2. Compute Kij for i ^ j. Using the results from Step 1 in equation (IX.2.6) 
gives ' 


N 


K ‘i = Yf, 'EiW” + # j. 


k = 1 


(IX. 2. 19) 


Step 3. Determine the eigenvalues and eigenvectors of the matrix K in equation 
(IX.2.10). H 

Step 4. Evaluate 11 and then R using equations (IX.2.14) through (IX.2.16). Let R n 
and K n denote the R and K matrices obtained in the n th iteration. K n is updated by 
using equation (IX.2.9) where R is replaced with R n . 
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Step 5. Check convergence. If 


||R„ - R n _i|| < S R , or ||K„ - K n _!|| < S K , (IX.2.20) 

then stop the iteration, otherwise, repeat the iteration from Step 3. In equation (IX.2.20) 

i? aiid Sr are two given small numbers, and || • j| denotes the norm of a matrix. The 
infinite norm 

IIRIloo = max |R 0 | (IX. 2. 21) 

is used in this study. 

It should be pointed out that in additional to matrices R and K, the eigenmatrix 
$ and eigenvalues yf are also obtained at the end of the above iteration. These are the 
results we need for the subsequent fatigue analysis. Most applications of the method of 
equivalent linearization end at yielding the matrices R and K. 

We have thus far completed the description of the method of equivalent linearization 
or determining the response of the nonlinear beam under consideration. The general idea 
of this method applies equally well to other nonlinear structures. To analyze the fatigue 
damage of the beam, we need samples of the time history of «,•(*) to generate time histories 
o the strain or stress. It turns out that it is more efficient to do so if the strain is expressed 
in terms of modal coordinates /?,(*). Using equations (IX.1.3) and (IX.1.4), the maximum 
strain of the beam at x = 0 or / can be written as 



(IX. 2. 22) 


where 


a = 



h &U 0) 1 

2 dx 2 /’ 


(IX.2.23) 


and I is defined in equation (IX.2.7). By substituting the first of equations (VI.12) into 
equation (IX.2.22) the maximum strain may be expressed in terms of the /?,(<), 


e = a T l + (IX. 2. 24) 

where 

a = a T $, f=$ T I$. (IX. 2. 25) 

The modal responses of the equivalent linear system, /?,(<), can be readily simulated 
from equation (IX.2.13). Equation (IX.2.24) then gives the strain time history. Once the 
strain time histones are generated, a fatigue damage model may be applied to estimate 
the fatigue life. 

efficient way of generating the time history of is to integrate equation 

(IX. 2. 13) numerically using the central difference scheme: 


W»+0 = Ci/3,'(<n) + C 2 A(f n -i) + C 3 fiW 0 (t n ). (IX. 2. 26) 
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where t n is the value of t at the n th time step, t n = n- Ai, where At is the time increment. 
The constants Ci, C2, and C3 are given by 


Ci 


2 - 7, 2 A t 2 _ - 1 ca<2 _ j 

1 > '-'2 — ■j , O3 — -7 . 

2 <A# + 1 -(A< + 1. 2 CA<+1 


(IX. 2. 27) 


Note that because the calculations are performed in terms of uncoupled modal coordinates, 
the modal responses may be obtained using parallel processors with vectorization. 

The local truncation error of the central difference scheme is of the order 0(At 2 ). The 
spectral stability condition of the scheme for each mode is given by [l] 


0< A<< 2-, i = 1,2,. ..,7V. (IX.2.28) 

Hence, the stability condition for the entire system is given by 

0<Af< A (IX.2.29) 

where y N is the largest eigenvalue in equation (IX.2.11). The time step A t in our simulation 
is chosen according to equation (IX.2.29). The time series of the white noise excitation 
W 0 (t n ) is simulated by using Rice’s approach and the Fast Fourier Transform [2]. 

We shall refer to the method of fatigue analysis described in this section as the quasi- 
simulation approach. In summary, the method of equivalent linearization is used to obtain 
the eigenvalues and eigenmatrix of the uncoupled equivalent linear system, and then the 
time histories of the modal coordinates are simulated numerically to generate the strain 
time histories for the fatigue analysis. It should be noted that once the equivalent lin- 
ear system is determined, the computational effort of simulating the time series of the 
response is identical to that required for a linear system and is independent of the type of 
nonlinearity. Since a large number of very long time series are required to obtain accurate 
fatigue estimates, the over-head’ required in constructing the equivalent linear system is 
comparatively small. 


IX. 4 Numerical results 

The quasi-simulation approach developed in the previous section has been applied to 
estimate the fatigue life of a nonlinear beam and the results compared to those obtained 
using a conventional numerical simulation. The comparisons were performed to identify 
the errors in the approximate method and to quantify the savings in computation time. 
The conventional simulations were obtained using a fourth order Runge-Kutta algorithm. 
The same time step, At, was used in all computations. The parameters of the beam used 
in the numerical study are: 

E = 10 7 #/m 2 , p = .l#/m 3 , h - .032 in, 

1 = 15m, T) = 3.016, b=2, A = .064, 
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where the variables are defined in the discussion following equation (IX. 1.1). The number 
of resonant modes was taken to be three. 

Estimated fatigue lives are shown in figures IX.4.1 and IX.4.2 for the nonlinear beam 
modeled in the previous sections. In these figures the dots indicate results obtained using 
a full numerical simulation and the lines are the result of the approximate method. The 
estimated fatigue lives are found to be in excellent agreement for a range of excitation 
levels and for a range of values of the fatigue exponent, 6, in equation (111.1,1). 

In the results shown in figure IX.4.1, damage was assumed to occur at each stress 
maximum in the time domain. In figure IX.4.2 damage was accumulated only for the 
greatest stress maximum between zero crossings. Of course, as discussed above, one can 
accommodate any desired damage counting scheme in a simulation procedure. 

Equation (III. 1.1) was implemented in the simulations in the form 

ISI 6 

D m «— D m + (IX.4.1) 

where D m is initially set to zero and S is the level of the stress maximum. The absolute 
value of the stress is taken to allow for positive contributions to damage due to negative 
maxima, b is the fatigue exponent which for most materials varies between 2 and 6. c is a 
material constant taken to be 6.56 x 10 30 . The value of c has no effect on the comparison 
between the methods shown in figures IX.4.1 and IX.4.2. 

Once the simulated damage is accumulated according to equation (IX.4.1) for a suffi- 
ciently long time, r, the average damage rate is simply 

A = D m /r, (IX.4.2) 

and the simulated mean fatigue life is 

T = 1/A. (IX. 4.3) 

Figures IX.4.1 does not reveal the amount of nonlinearity in the system at the excita- 
tion levels used here. To indicate the influence of nonlinearity, the power spectral density 
of the stress time history was computed for three excitation levels. The results are shown 
in figures IX.4.3 through IX.4.5. In figure IX.4.3 it is clear that three resonant modes 
contribute to the response. In this figure the excitation level was not intense enough to 
elicit nonlinear behavior. As the excitation level is increased, the full simulation shows 
that the three dominant peaks are shifted to higher frequencies and become broadened. 
Additional peaks in the stress spectrum also become apparent at the higher excitation 
levels. Since the response is drastically different at the higher levels than at low excitation 
levels, nonlinear effects clearly play an important role. A further indication of the degree 
of nonlinearity in the response can be seen in figure IX.4.6 which shows the mean square 
displacement at the center of the beam with and without nonlinearities included. The 
figure shows that above roughly -35 dB (g 2 /Hz) the nonlinearities are substantial. 

The increase in resonant frequencies and the broadening of the resonant peaks are 
common behaviors of random systems with stiffness nonlinearities. The broadening of the 
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resonant peaks has been shown to result from the fact that the resonant frequencies become 
random in a system where the stiffness depends on the random response amplitude [31. It 
has also been shown (and can be seen in figures IX. 4.3 through IX.4.5) that while equivalent 
linearization accounts for the increase in the resonant frequencies due to nonlinear stiffness, 
it does not depict the broadening behavior [4]. Although equivalent linearization does not 
accurately represent the response spectrum, it has been applied very successfully using a 
probabilistic approach to estimate fatigue life for a single degree of freedom system [5], 
The results obtained in the present study show that equivalent linearization may also 
be used in a time domain approach to estimate fatigue lives for multi-degree of freedom 
nonlinear systems. Theoretical aspects, such as error analysis and convergence, of the 

method of equivalent linearization for the fatigue analysis of nonlinear structures await 
further studies. 

Along with the shift and broadening of the resonant response peaks, figures IX.4.3 
through IX.4.5 also show that the nonlinearities produce additional peaks in the stress 
spectrum. These added peaks arise from the nonlinear stress/displacement relationship in 
equation (IX. 1.2) and from the nonlinear coupling in equation (IX. 1.6). 

As mentioned above, the numerical comparisons were performed to evaluate the accu- 
racy of the approximate simulation and to determine the improvement in computational 
efficiency relative to conventional simulation methods. In the present effort, our purpose 
is not to conduct a comprehensive comparison of numerical simulation methods The 
present results primarily serve to indicate the validity and potential of the approach. For 
the cases studied here where the beam was assumed to have three resonant modes the 

quasi-simulations required roughly (l/8)<* the computer time used to perform the full 
simulations. 

All of the calculations were performed on an IBM 3090 computer using a single proces- 
sor without vectonzation. Since the present approach is well suited to parallel computation 
on multiple processor computers, even greater savings are achievable. 
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Figure IX.4.4 

Comparison of stress spectra obtained using a full 
numerical simulation and the present quasi— simulation approach, 
g is the acceleration due to gravity. 
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X. Nonlinear Plates 


In the previous section an application of equivalent linearization is developed to esti- 
mate the fatigue life of a simple beam with multiple resonant modes. The results indicate 
that the method is applicable to fatigue prediction in nonlinear structures. Because we 
are usually concerned with the problem of estimating the fatigue life of structures that are 
considerably more complicated than beams, the procedure has been applied to a nonlinear 
plate. The next logical application of the method should be on stiffened plates or shells 
such as fuselage structures. This extension of the method will be left to future studies. 

In the present section, a model of a nonlinear plate is developed based on Von-Karman 
plate theory and a governing equation for transverse motion only is derived by neglecting 
the in-plane inertia. As in the analysis of the nonlinear beam of the previous section, 
an equivalent linearization approach is used to approximate the nonlinear system with a 
equivalent linear system by selecting the damping properties and natural frequencies to 
minimize the mean square error between the approximate and original system. Further, the 
construction of the time domain response due to stationary Gaussian random excitation is 
relatively straight-forward and the geometrical nonlinearity which results in a nonlinear re- 
lationship between strain and displacement is considered. Because the plate we considered 
has multiple resonant modes, the response will be broad-band in nature. The Rain-Flow 
Cycle Counting scheme which is applicable to a broad-band process is used to estimated 
the fatigue life. The present approach is based on an analysis of the stress and strain in 
the time domain because the Rain-Flow Cycle Counting scheme is extremely difficult to 
incorporate in a probabilistic formulation. 

It has been shown previously that although the equivalent linearization method is not 
appropriate for estimating the power spectral density and some other statistics, it yields 
reasonably accurate estimates of the mean square response. We have shown that the 
estimated fatigue lives of nonlinear beams are in excellent agreement with results obtained 

by more accurate techniques. The purpose of the present study is to extend these results 
to nonlinear plates. 

Comparisons are presented of the estimated fatigue lives of a plate with multiple 
resonant modes and geometrical nonlinearity as obtained by present method and by a 
full numerical simulation. It is found that the fatigue lives agree very closely even when 
nonhneanties dominate the response. Since the major computational effort of the present 
approach is concentrated on constructing the time domain response of the equivalent linear 
system, the efficiency and complexity of the analysis are not affected by the nonlinearity 
of the structure. This feature makes the present approach a very practical method for 
fatigue study of highly complex structures. Further, the equivalent linearization approach 
is efficient enough to allow us to compute fatigue lives at a large number of points on a 
plate which is almost impossible to get by using full numerical simulation. 

A comparison is presented of predicted fatigue lives with and without nonlinearities 
included in the model. It is found that although the predicted fatigue lives are in close 
agreement when excitation is low, the fatigue lives differ significantly when the excitation 
is high. This is because the geometrical nonlinearities which are caused by the stretching 
of the plate reduce the mean square stress response and increase the predicted fatigue 
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lives. It is implied that we can not calculate the low excitation fatigue lives from the high 
excitation fatigue lives by using a simple linear extrapolation. 

X.l. Nonlinear Plate 

In this section, the governing equation for a nonlinear plate which is subject to random 
excitation is obtained, an equivalent linear approach is derived in detail and the maximum 
stress of the plate is computed in terms of the modal coordinates. It will be assumed that 
the plate is excited by the transverse motion of its boundaries. 


X.1.1 Basic Assumptions for a Nonlinear Plate 

The assumptions of the present investigation are based on the Von-Karman plate 
theory. For a thin plate subject to large elastic deflection with small rotation, the following 
assumptions are made: 

1. The plate is thin, the thickness h is much smaller than the wavelength of any vibration 
modes considered in the model. 

2. The magnitude of the deflection W is of the same order as the thickness h of the plate. 

3. The slope is everywhere small. 

4. The in-plane displacements U and V are infinitesimal. In the strain- displacement re- 
lations, only those nonlinear terms involving are included, all other nonlinear 

terms are neglected. 

5. All strain components are small so that Hook’s law holds. 

6. Kirchhoffe s hypotheses hold: Tractions on surfaces parallel to the middle surface are 
negligible, and strains vary linearly within the plate thickness. 

7. The in-plane inertia is negligible compared with the transverse inertia. 


X.l. 2 Governing Equation for a Nonlinear Plate 


Let U(x, y,t), V(x, y, <),and W{x, y, t)(see Figure X.l) be the displacements of a point 
(x,y) at the middle plane of the plate relative to the motion at the boundaries at x = 
0, a, y = 0, b. The displacements of an arbitrary point within the plate are given by, 


U x (x, y, z,t) = U(x,y,t) — 


U y (x,y,z,t) = V(x,y,t) - 
U z (x,y,z,t) = W(x,y,t). 


d Wjc^t) 

Z dx 

dWQ,y,<) 

Z dy 


The Green’s strain components are given by, 


* _ j i dU j , x^dUkdUk^ 
” r dx, + dz, ^ Ox, dx, 1 ' 


(X.1.1) 


(X.1.2) 
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where we have dropped the explicit dependence on x,y and t to simplify the notation. 

Because we regard the thin plate as a plane stress problem, we can neglect the z 
components of the strain and stress. From equation (X.1.2) and Hook’s law, we can get 
the relationships for strain e and stress a in terms of the displacements at the middle plane 
of the plate, 


£xx 


e yy — 


dU_ _ d 2 W 1 fdW\ 2 
dx Z dx 2 + 2\lh) 

dV_ (FW 1 / dW \ 2 
dy Z dy 2 + 2 \ dy ) 

1 ' 


e xy — - 


du dV <FW dW dW 

2 [ dy dx Z dxdy ^ dx dy J 


(X.1.3) 


0TT. = 


E 


1-1/2 

E 


( C XZ + Vtyy) 


a yy ~ i ~ ( e vy + vt *x) 


l-I/ 2 

E 


a *y ~ 


l + i/ 


(X.1.4) 


where E , v are Young s modulus and Poisson ratio of the material, respectively. 
The kinetic energy E and strain energy V are given by, 


h. , 

2 6 a 


-?///(£ )♦(*)♦(£ 

k r\ a ' ' ' 


+ 


-i 0 0 

T b a 


V = 


-!///< 

0 0 

2 


^yy^yy H” xy^xy} dxdydz , 


dW 0 \ 2 ] , , 


Ot J dxdydz 

(X.1.5) 

iz, 

(X.1.6) 


where, a, 6, /i are the length, width and thickness of the plate, p is the mass density and 
W 0 is the prescribed motion of the boundaries. Substituting equations (X.1.1), (X.1.3) 
and (X.1.4) into (X.1.5) and (X.1.6), integrating with respect to z, gives, 


-+/[( 


dU 

dt 


o o 

6 a 


avy /gw sw„\ 2 i 

~m) + {w + ~dF) \ dxdy 


+ (an 


+ 


ph 3 f [\(d 2 W\ 2 /d 2 W\ 2 ] , 

W/ [VaJsj + [Wt) j y; 


0 0 


(X.1.7) 
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b a 


5° 

ii 

(d 2 W 
A 9x2 

\ 2 fd 2 W' 

J Iv, 

n d 2 Wd 2 W 
) +2t/ dx 2 dy 2 

+ 2(1 - v) 

(d 2 w\ 2 ) 

\dxdy) 

dxdy 

6 a 

♦?. // 
0 0 

(i) 

'•©'-SI 





1 - 1/ 

/fdUV , 

fdV\ 2 n dUdV 

w 



+ 2 

VlaA + l 

[dxj + 2 dy dx 

J dxdy 


(X.1.8) 

b a 

*!//[ 
0 0 

fdW\ 

\dx) 

'*(?)' 

4.9 f<W\ 2 fdW 

\dx J \dy 

^ J dxdy 




+ 


- t f\ dU ( dW \ 2 I 9V(dW\ 2 \dU ( dW\ 2 dV(dW' 

^ J J [d x \ dx J dy \ dy J [dx \ dy ) ^ dy\ dx / 


2i 


Jdu x dv\awdw 

+ ( 1 “ I/ ) -5T + -5T 5 - 


V dy dx J dx dy 


dxdy, 


where, 


D = 


Eh 3 

12(1-1/2) 



(X.1.9) 


The second term of (X. 1.7) represents the rotational inertia and it will be neglected 
because it only influences the very high frequency responses. This will cause the mass 
matrix to be diagonal. The first term of (X.1.8) will generate the linear term of the 
stiffness matrix of the transverse motion, the second term will generate the linear term 
of the stiffness matrix of the in-plane motion, the third term will generate the nonlinear 
term of the stiffness matrix of the transverse motion and the fourth term accounts for the 
coupling between transverse motion and in-plane motion. 

The virtual work done by the non- conservative forces may be written as, 


b a 

SW = -J ' J c( Xt y)?Kew dxdy, 
0 0 


(X.1.10) 


where c(x, y) is the viscous damping coefficient. 

Now, we assume the displacement responses U, VkW can be expanded in a finite set 
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orthogonal functions, 


M 

U(x,y,t ) = ^2ai(x)f3j(y)Uj(t) 

i=i 

M 

V(x,y,t) = 5^7 i(x)nj(y)Vj(t) (X.1.11) 

7=1 

N 

W(x,y,t ) = ^2 <f>i(x)ip j (y)W I (t), 

1=1 

where a, (a:), /?y(y), 'Ji(x), r)j(y) are the functions used to describe the in-plane motion, 
4>i{x), ipj(y) are the functions for the transverse motion, M, N are the number of terms 
used for the in-plane and transverse motion, respectively, and ~ denotes the index relative 
to the in-plane motion. Uj(t), Vj(t), and Wj(t) are unknown functions of time. i,j are 
selected so that any desired number of half waves is included in each direction. The 
relationships between the indices are 


and, 


I = (i - 1) x rij + j 

for i = 1,2, j = 1,2, 

I = (i - 1) x hj + j 

for i 1, 2, • * • n j, j — 1, 2, • • • 7 i j 


(X.1.12) 


where rii,rij are the number of modes for transverse motion in the x and y direction, 
respectively, and hi, rij axe the number of modes for in-plane motion in x and y direction, 
respectively, and 


N = 71,' X fly 

M = hi x hj. 


(X.1.13) 


An example is given in table X.1.1 
with hi = hj = 3, 

i i 


1 

2 

3 

4 

5 

6 

7 

8 

9 


1 

1 

1 

2 

2 

2 

3 

3 

3 


to show the relationship between I and i,j for U(x,y,t ) 

j 

1 

2 

3 

1 

2 

3 

1 

2 

3 


Table X.1.1 An Example of Index Relationship 

In the present study, the boundary conditions for in-plane motion are taken to be 
fixed and for transverse motion they are assumed to be clamped-clamped. A reasonable 
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choice for the shape functions used in equations (X.1.11) is the set of eigenfunctions for a 
fully fixed beam, 

, S . 2l7TX 

Qri( x J = 5m ' 


a 

jwy 


Pj{y) = sin 


e x . 17TX 

7 i(xJ = sm 

a 


r h(y) = sin 


2j*y 


— cos( —■) - cosh(^—) + Di(sin (- - smA(— )) 

<h(v) = cos ( - cosh(^-) + _ 5 i„A(M)) j 

where p, and Z), axe given in the table X.1.2 for the first six eigenfunctions. 


(X.1.14) 


2 

1 

2 

3 

4 

5 

6 


Pi 

4.730040745 

7.853204624 

10.99560784 

14.13716549 

17.27875996 

20.42035225 


Di 

-0.982502215 

-1.000777312 

-0.999966450 

-1.000001450 

-1.000000000 

-1.000000000 


Table X.1.2 Eigenfunction Coefficients for a Clamped- Clamped Beam 
Introduce following notations, 


b a 

M(I,J ) = ph J J dxdy 


o o 

b a 


<500 — phj j dxdy 
o o 

b a 

mu(I,J ) = ph j J (atiftjockfli) dxdy 


(X.1.15) 


o o 

b a 


m v (i,J) = phj J (jir)jy k r)i) dxdy, 
o o 


b a 

k w (I,J) = dJ JtfiWfri + WiMi + 2^^) dxdy 
0 0 
b a 

ku{I, J) =■ B J J — — atiPjOikPi) dxdy 

o o 
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b a 

kv(i, J) = B J J (yiVjlkrj'i + ~2~7i r )j7k r )i) dxd V 
0 0 

b a 

k uv (i,j) = J Zva'iPjKkVi + (1 - u)a iP'j~i' k m) dxd V (X.1.16) 

0 0 

b a 

k WW (I,J,K,L) = ^ J J (<f> , i Tpj4>' k *{’l<t>' m rl>nK4 , s + </>i^j<i>k^l<f>m‘tp' n <f>ri’s 

0 0 

+ 2<t>' i xl>j<f> , k ipi<f> m xj>' n $ r ij;' a ) dxdy 

b a 

kwu(K,I, J) = B j J (tiipj^tpta'nPn + vW j M , l a' m P n + (1 - dxdy 

0 0 
b a 

k WV (K,I, J) = B J J (WjMllf mTj'n + Vp\^j<i> k ^ llm Ti n + (1 - v)^ i4>k'l’Wm r }n) dxd V , 

0 0 


C(I, J)=j J c(x,y)</>iipj<f> k rpidxdy, (X.1.17) 

o o 

where f denotes the derivative with respect to the spatial variable. Then the energy 
expression can be written in a compact way. Substituting equations (X.1.11) through 
(X.1.17) into equations (X.1.7) and (X.1.8) gives, 

. N N 

T = o ( X X M ( ] ’ 

7=1 7 = 1 
M M 

+ E 3)Vfij + rn v (i, J)VfVj) (X.1.18) 

7=1 7=1 
N 

+2^<5(7)W 7 Wo + abW*) 

7=1 

1 N N 

'’^(EEww 

7=1 7=1 
M M 

+ E E(M J ' WjUj + k v (i, J)V,V } + 2 k uv (i, J)UjVj) 

7=1 7=1 

- N N N N 

+ 2^EE £ W(J, l K, LWjWjWkWl (X.1.19) 

1 1=1 J= 1 K= 1 L= 1 
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M N N 

+ 2 T,T,^ k »'v(K,I,J)W I W J U lt: + k W v(k,I,J)W I WjV it )), 

k= i *=i J=i 


where we have dropped the explicit time dependence of t/,j(<), and W, } (t) to simplify 

the notation. 

The Ritz Method will be used in this model to find the mass and stiffness matrices 
and the response of the plate by using Hamilton’s Principle, 

*2 

6 J[T - V + W] dt = 0, (X.1.20) 

<i 

where 8 is the variational operator. Note that the base excitation, W 0 , is a known function, 
so that, 

<5Wo = W 0 = 0. (X.1.21) 

Applying the variational operation to equation (X.1.18) and (X.1.19) gives, 

H M 

^ = (53 + Q(I)W 0 )8Wj + J)Uj8Uj + m v (I, J)Vj8(ffi, 1 . 22 ) 

J=1 7=i 

N N N N 

6V = ( 53 M A J)Wj + £ 53 53 k ww(I, j, K , L)WjW k Wl 

J= i j=i jf=i l=i 

N M 

+ 2 53 53 {kwv{K , I , J)WjUfr + kwv(K,I, J)WjVj{))8Wij 
J=1 k= 1 

M M N N 

+ (53 Wj + 53 £ /)v> + 53 53 K ) w jw k )^h,i.23) 

J=1 7=1 J=1 A"-l 

M M N N 

+ (53 m £ 'W + 53 wc/, + 53 53 

J=i j=i j=i 

Substituting equations (X.1.11) and (X.1.17) into equation (X.1.10) gives, 

N 

6W = - 53 C(I, J)Wj8W T. (X.1.24) 

J= 1 


Substituting equations (X.1.22) through (X.1.24) into (X.1.20), we can get the governing 
equations for Wj, and Vj in the expansions in equations (X.1.11), 

53 m ( j , + 53 C ( J ’ + 53 kw ( J > *0^ + 53 53 53 *ww(/, j , k, l)Wjw k w l 

J=1 J= 1 J=1 J=l*T=lL=l 
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(X.1.25) 


N M 


+ 2 Y^Y,( kw v{KJ,J)U i Wj + k wv (KJ,J)VKWj) = -Q(I)W 0 , 

J=1 K=1 

for I = 1,2, • •• , AT 


M 


M 


M 


7= 1 


M 


7= 1 


7=1 


M 


7=1 


7=1 


M 


7=1 


N N 

^m(/, 7)£/j + ^ &[/(/, J)C7 j + ^ kuv(I, J)Vj + ^ ^ k wu (I, 7, K)WjW k (X.Q,26) 

7=1 7T=1 

for 7 = 1,2,---,M 

^m(J , J)Vj + ^ &v(7, 7 )Vj + ^ kuv(I, 7)Vj -f ^ ^ kwv{I, 7 , TQlFjrW/f 4^CL.27) 

for 7 = 1, 2, 


AT IV 


7=1 A'=l 


From previous analysis, we know [M(7, 7)] is a diagonal matrix having diagonal ele- 
ments M/, if we assume c{x , y) is a constant, [C(7, 7 )] will also be a diagonal matrix having 
diagonal elements C/. In the present study, the general terms of the stiffness matrix have 
been evaluated by using the symbolic computer program MACSYMA. 

According to the assumption of the Von-Karman plate, we can neglect the in-plane 
inertia, UjandVj and solve for Z7j, Vj by using (X.1.26),(X.1.27), 



[£[7(7)7)] 

[kuv( 7, 7)] 


[kuv{1, 7)] 
(M7,7)] 


( 


N N 


{ E E 

K= 1 L= 1 




N N 


(E E 

K= 1 L = 1 


kwu(J, K, L)W k Wl] ^ 
kwv(j,K,L)W K W L ] 


J 


(X.1.28) 

where {17(7)} is a vector having elements 17(7) for 7 = 1,2, • • • , A7, etc. and [ku(1, 7 )] is 
a matrix having elements kjj( 7, 7) 


Equations (X.1.25) through (X.1.28) may be manipulated to obtain the governing 
equation which depends on the transverse displacement W only. 


N N N N 

M/W/ + + J2k N (I,J,K,L)WjW K W L = -QiW 0 , 

7=1 7=1 K = 1 L= 1 

(X.1.29) 

where kw{1 , 7) is the linear part of the stiffness matrix and £jv(7, 7, K , L) is the nonlinear 
part. 


X.1.3 Equivalent Linearization 

In order to predict the fatigue life, a sufficiently long time history of the response under 
random excitation is needed. Because of the nonlinear term in (X.1.29), it is not practical 
to solve for the response in the time domain using full numerical simulation methods if 
N is large. The numerical effort in evaluating equation (X.1.29) is roughly proportional 
to AT 4 . In this section, we shall present the method of equivalent linearization for solving 
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equation(X.1.29). This will lead to an approximate solution to equation (X.1.29) in which 
the numerical effort is proportional to N 2 rather than N -4 . This makes it practical to 
include a large number of terms. 

The method of equivalent linearization consists of approximating equation (X.1.29) 
by an equivalent linear system, 


1 N 

w. I + tiWj + — Y W J)Wj + ei(W) = q(I)Wo , (X.1.30) 


j=i 


with 


f = £l 

^ M t 




Qi_ 

Mi 


(X.1.31) 


where k e (I,J ) is the equivalent linear stiffness matrix and ej(W) is the error term which 
can be written in a vector form, 


e = {ej(]£)} 

N 


f N N N N N 

= 1 W J)Wj - k w(I, J)Wj -YYsYj M J > J, K, L)WjW k W l 1, 

kj=1 J=\ J=1K=1L=1 j 


(X.1.32) 

where the index I indicates the row of the vector. 

The elements k e (I, ,7) are chosen to minimize the steady-state mean square value of 

e, 


3E [e^e] 
dh(I, J) U ’ 

where E[-] denotes the expected value. 

Expanding (X.1.33), and using (X.1.32) give, 


(X.1.33) 


N n 

Y K(I,M)E[W m Wj]= Y kw(I,M)E[W M Wj] 

W=1 M=l 

N N N (X.1.34) 

M= 1 K=1 1=1 


for I = 1 , 2 , • • • , N. 

Now we assume Wo is Gaussian white noise, so the solution for the Wi is known 
to be jointly Gaussian and all the odd order steady-state moments of the Wj are zero. 
The fourth order moments in equation (X.1.33) can be expressed in terms of lower order 
moments, 


E[WjW m W k W l ] = E[WjW m }E[W k W l }+ 
E[WjW k )E[W m W l ) + E[WjW l ]E[W m W k ]. 
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Substituting (X.1.35) into (X.1.34) gives, 

N N 

k c (I,J) = k w (I,J)+Y, ^ kN (^ J ^ R ^ L ) + 2k N(^K,L,J))Rw(K,L), (X.1.36) 

K= 1 L = 1 


where, 

Rw(K,L ) = E [W K W L ] (X.1.37) 

is the steady-state correlation matrix for Wjc, K = 1, 2, • • • , N. 

Equation (X.1.36) can be used to determine the components of (/,/). However, 
one has to solve for [fc e ] and [f?vv] simultaneously. An iterative procedure is described in 
the following. Let $(/, J ) be a matrix whose columns are eigenvectors for 



(X.1.38) 


«^) is a normalized modal column matrix with the following properties, 


where, 


1 = 


/I 0 

0 1 

Vo o 


[*F[*] = I 
PffiAf]- 1 [*,][«] = a 2 , 


°\ 

0 


1 ) 


n 2 = 


(u\ o ... o 
. 0 


0 U>2 


Vo o 


UJ 


N 1 


(X.1.39) 


(X.1.40) 


and uij is the Jth natural frequency. Define a transformation to modal coordinates, Bj, 


N 

Wi = Hr, J)Bj 

V (X.1.41) 

// = £ Hr, J)U- 

J= 1 


Substituting (X.1.40) into (X.1.30) and neglecting the error term leads to an uncoupled 
set of differential equations, 


Bi + tiBi +u}Bi = fiW 0 . (X.1.42) 

Assume Wo{t) is Gaussian white noise with single-side power spectral density G^. We 
define the steady-state correlation matrix for the modal coordinates J3/, 

Rb(I,J) = E[B!Bj]. (X.1.43) 
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From equations (X.1.41) and (X.1.42) it can be shown that 


E[BjBj] 


UifjGw, 

("? + 2 ^+^)' 


(X.1.44) 


The correlation matrices [.Rvy] and [f?#] having elements -Rw(/, J) and _Rg(/, J), respec- 
tively, are related by 


[Rw] = [$][fl B ][$] T . 


(X.1.45) 


Using equations (X.1.38) through (X.1.45), the iterative procedure to determine [jfc e ] 
and [iliy] may be stated as follows: 

Step 1. (Zero order approximation) Assume that 


Rb{I, J ) = 0. 


Then 


R\v(I , J) = 0 
k e {I,J) =kw(I,J)- 


(X.1.46) 


(X.1.47) 


Step 2. Determine the eigenvalue and eigen matrix of equation (X.1.38). 

Step 3. Evaluate R B and then R w using equation (X.1.43) through (X.1.45). Then [jfc e ] is 
updated by using (X.1.36). 

Step 4. Check convergence. We will stop the iteration if the convergence conditions are satis- 
fied, otherwise, repeat the iteration from Step 2. 

It should be pointed out that in addition to the matrices [fZw]» [frs] and [jfe e ], the 
eigen matrix $ and equivalent linear natural frequencies u > / are also obtained at the end 
of the above iteration. These are the results we need for the subsequent analysis. 


X.1.4 Stress Calculation 


In order to predict the fatigue life, we need the stress or strain expressions in terms of 
the equivalent linear modal coordinates. The stress and strain should include both linear 
bending and nonlinear stretching components. Substituting equation (X.1.3) and (X.1.11) 
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into equation (X.1.4) gives, 


E r ~ - 

* 1=1 

N 

~ z + v<i>i( x )i’"(y))Wi 

7=1 

^ N N 

Z 7=1 7=1 
£ M 

^ = TZr^^^ x ^W v l + va'iWPjiypj) 

X ^ 7=1 

77 

~ z (y) + rf(aO^(y))W7 

7=1 

j iV iV 

+5 ED* ,(x)V-' (y)Mx)rp't(y) + ^(^)^(y)^(^)V’i(y))W / Wj) 

7=1 7=1 

E M 

a *y = 2(1 + JD Q2^ x ^(y) u i + 7i( x ) r h(y) v l) 

N 

-2z^$(x)V>'(y)W7 

7=1 
77 77 

+EE ^.(x)^(y)^ fc (x)V>;(y)W7W7). 

7=1 7=1 

(X.1.48) 

From equation (X.1.48), we know that the stress of the plate is multiaxial and an 
equivalent stress approach should be developed to calculate the fatigue life. But, because 
the equivalent stress has different directions at different points, it is difficult to construct 
and count the stress cycles when using a cycle counting scheme to predict fatigue life. 
According to the basic plate theory, the maximum stress will occur at the boundary if the 
boundary condition is clamped and either a xx or a yy will dominate along the boundary. 
So in this project we simply calculate the fatigue life by using a xx and a yy separately. 
Along the boundary, the maximum a xx will occur at the middle of the edge in y direction 
while the maximum a yy will occur at the middle of the edge in x direction. If we ass um e 
the width b is greater than the length a, then 


&max — & xx (x — 0,y — — ,2 — 


(X.1.49) 


where ± means we can choose either the upper plane or bottom plane of the plate and it 
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will come up with the same solution. Define 


SU(!) = W><\) 

S v{I) = 

S w( T ) = 2 ( 1 ^) (^ 4 ) + 

Sww(I,J ) = 

Substituting equation (X.1.48) and (X.1.50) into (X.1.49) gives, 


M N N N 

Vmax = ^2{Su(I)Uj + fiSy(i)Vj) ± ^ Sw(I)Wi + ^ ^ Sww{I , J)WjWj. 

1=1 1= 1 7=1 7=1 

(X.1.51) 

Recalling the modal transformation equation (X.1.41), we can express equation (X.1.51) 
as, 

M N N 

a max — / + flSv(i)Vj) ± <SW(J)$(J, J)B J 

1=1 1=1 7=1 (X.1.52) 


N N N N 
7=1 7=1 K= 1 L=1 


where Vj are given by equation (X.1.27) and can also be expanded as an expression in 
terms of the equivalent linear modal coordinates £/. 


X.2. Quasi-Simulation and Fatigue Estimates for Multi-mode Nonlinear Plate 

As in previous sections, we will assume that the relation between the stress amplitude, 
5, and the number of the cycles to failure, N, is given by the well known equation, 

N = c/S\ (X.2.1) 

where c and b are experimentally obtained constants for a given materials. For most 
material b vanes between 2 to 6 and c is taken to be 6.56 x 10 30 here. For the purpose of 
estimating high cycle random fatigue life, the Palmgren-Miner linear damage accumulation 
rule may be applied. In this theory, the total damage D, is written as the sum of the damage 
due to all damage events, 

D = J2 AD i- (X.2.2) 

i 

Failure is predicted to occur when D = 1. The fatigue lifetime is then the amount of time 
it takes for this to happen. 
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The damage increment due to one cycle, AD, is given by, 

AD = \Si\ h /c, (X.2.3) 

where S{ is the stress amplitude of the ith damaging event determined by the cycle counting 
scheme. In our study, because the band width of the response of the multi-modes nonlinear 
plate is very wide, we use Ram-flow Cycle Counting scheme. 

Equation (X.2.2) may be implemented in a simulation of the damage accumulation in 
the form 

D <— D + (X.1.4) 

where D is initially set to zero. Once the simulated damage is accumulated according to 
equation (X.2.4) for a sufficiently long time, r, the average damage rate is 

A = D/t, (X.2.5) 

and the simulated mean fatigue life is 

T = 1/A. (X.2.6) 

Before counting the damage, a sufficiently long applied force time history should be 
simulated, then the time history of the modal responses Bj are calculated using equation 
(X.1.42) by the central difference method. Finally, the time history of the stress is calcu- 
lated by equation (X.1.52). It should be noted that once the equivalent linear system is 
determined, the computational effort of simulating the time series of the response is iden- 
tical to that required for a linear system, and is independent of the type of non-linearity. 
So the computational time is significantly reduced compared to the direct simulation of 
the nonlinear system. 


X.3. Numerical Results 

The approach developed in the previous section has been applied to estimate the 
fatigue life of a nonlinear plate, and the results compared to those obtained using a con- 
ventional numerical simulation. The comparisons were performed to identify the errors in 
the approximate method and to quantify the savings in computation time. The same time 

step, At, was used in all computations. The parameters of the plate used in the numerical 
study are: 

E = 10 7 #/in 2 , p = 0.1 #/m 3 , h = 0.032 in, 

a = 10 in, b = 15 in, p = 0.3, 

where the geometry of the plate is shown in Figure X.l. The number of resonant modes 
for transverse motion, N , was taken to be 9 and the number of resonant modes for in-plane 
motion, M, was taken to be 256. 
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The estimated fatigue lives for each point of the plate under different excitation levels 
are shown in Figures X.2 to X.4, which are calculated assuming that cr xx is the only 
significant stress component. The location of the minimum fatigue life is found to be 
at the middle of the longest edge, where the maximum stress occurs. Because of the 
computational time limitation, it is almost impossible to get a similar figure by using full 
simulation. The comparison of the approximate results and those obtained using the full 
system in equations (X.1.29) is performed at the point where the minimum fatigue life 
happens. 

Estimated fatigue lives at the middle point of the width under different excitation 
levels are shown in Figure X.5. Results obtained using the full simulation and the equivalent 
linearization methods are found to be in excellent agreement for a range of excitation levels. 
To indicate the influence of the non-linearity, results are compared in figure X.5 between 
a nonlinear plate and a plate which is assumed to respond linearly. The figure shows for 
input excitation levels Wo(t) above — 30dB g 2 /Hz, the plate response is highly non-linear 
and the mean square responses and fatigue lives deviate substantially from that of a linear 
system. This is because the geometry nonlinearities reduce the mean square values of the 
displacement and the stress and increase the predicted fatigue lives. So it is unreasonable 
to calculate the low excitation fatigue lives from high excitation data only by using a 
simple linear extrapolation. Another observation we found here is although we considered 
both bending stress and membrane stress in our stress model, the results from Figure X.5 
show that the bending stress is totally dominant at the middle point of the edges. We can 

understand it from equation (X.1.52) and recall that the boundary conditions we assumed 
are clamped-clamped. 

As mentioned above, along with the comparisons of the accuracy, the computation 
time was compared. In the present effort, the quasi-simulations required roughly one- 
hundredth of the computer time used to perform the full simulation. It is expected for 

more complex systems with more modes, the present method would prove to be even more 
advantageous. 

As expected, while the power spectral density of the stress has very good agreement 
between full simulation and equivalent linearization method at low excitation levels, they 
are totally different when excitation level is high. This is because we approximate a 
nonlinear system by a linear system and assume the response is Gaussian when excitation 
is. Gaussian white noise. Comparisons of the power spectrum of the stress are shown in 
Figure X.6 to Figure X.8. Although peak shifts are found and agree very closely from both 
methods, the peak broadness, which is a major feature of a high nonlinear system, is very 
obvious in the full simulation solution when excitation level is high but can not be found 
in the equivalent linearization solution. Also, some extra peaks are found in the equivalent 
linearization solution at high excitation levels at the points where resonant frequencies 
are doubled. This is caused by the nonlinear stress terms in equation (X.1.52). This 
phenomenon is not obvious in the full simulation solution because of the peak broadness. 
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XI.l Effect of Nonlinearity on Spectral Shape 

In previous sections we have examined the effect of spectral shape on the predicted fa- 
tigue lives of linear structures using several damage counting methods and have developed 
simulation methods to estimate the fatigue lives of nonlinear structures. It has also been 
pointed out that nonlinearities in a structure can have a pronounced influence on the re- 
sponse power spectrum in structures that are driven with random excitation. Although 
there is no simple way to estimate the fatigue life from the power spectrum of a nonlinear 
structure, it is important to understand the how the spectrum can be affected by nonlinear 
effects. In the present section, we present an approximate method of calculating the power 
spectrum of nonlinear systems. 

In experimental studies of the response of structures to random excitations it is very 
common to characterize the response by measuring the power spectral density. Vibration 
engineers are usually very familiar with this representation of the structural behavior and, 
in a glance, can gain considerable insight into the system being studied. One can eas- 
ily see how many resonant modes contribute to the response by looking at the resonant 
peaks in the spectrum and one can also determine whether the system is heavily or lightly 
damped. There are situations, however, where nonlinear effects in the structure influence 
the measured power spectrum in a manner that precludes straight-forward interpretation. 
In cases where the response has sufficient amplitude to elicit nonlinear behavior, the res- 
onant response peaks in the power spectrum can take a dramatically different form than 
that observed in linear systems. The main goal of the present investigation is to pro- 
pose a new approximate scheme for describing the influence of structural nonlinearities 
on the response power spectrum. It is hoped that the present study will provide a better 
understanding of the random response of nonlinear systems. 

The present approximate representation of the spectrum is applied to a bilinear oscil- 
lator in which the nonlinearity has a very pronounced influence on the response spectrum. 
Comparisons are presented of the estimated power spectral density obtained using the 
present approximate scheme and by direct numerical simulation of the governing equation. 
Excellent agreement is observed between the two methods. 

The effect of nonlinearities on the response power spectral density has been studied 
by a number of investigators. Early approximate analytical methods have been presented 
that are based on either a perturbation approach [1] or equivalent linearization [2,3]. These 
methods are found to give reasonable results only for very small nonlinearities. Other 
studies based on numerical simulations have shown that in the case of a Duffing oscillator, 
where the stiffness contains a cubic nonlinearity and the damping is assumed to be linear, 
the resonant response peak in the power spectral density tends to broaden and increase in 
frequency as the level of the random excitation is increased [4-6]. This behavior has been 
observed experimentally in studies of the high level random response of beams and plates 

[7,8]. 

An approximate analytical procedure was proposed by the author to estimate the 
response power spectrum of a Duffing oscillator with linear viscous damping [9], The 
method is based on an adaptation of the method of equivalent linearization where the 
resonant frequency of the equivalent linear system is assumed to be random. This pro- 
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vided very accurate estimates of the spectrum when compared to results obtained using 
numerical simulations. It has not been possible, however, to extend the method to more 
general nonlinear random systems. The main purpose of the present study is to propose 

an approximate scheme that may be applicable to a wider class of nonlinear systems than 
studied previously. 

The basic assumption of the method proposed in reference [9] is that the response 
spectrum is strongly influenced by nonlinearities in the system because in a nonlinear sys- 
tem, the stiffness, and hence, the natural frequency, depend on the response amplitude, 
n a system with sufficiently light damping, the random response will behave as a narrow- 
band process with an oscillation period that will vary randomly. The random variation 
of the oscillation period is a direct result of the fact that the response amplitude varies 
randomly. The random fluctuation of the dominant oscillation frequency will lead to a 
broadened resonant response peak in the power spectrum. It is reasonable to assume, 
therefore, that if certain statistics are known concerning the random amplitude fluctu- 
ations, and if the relationship between the system resonant frequency and amplitude is 
known, then one should be able to approximate the response spectrum. 

The procedure presented in the following is based on assuming that the nonlinear 
system may be approximated as a linear system having a natural frequency that varies with 
e response amplitude in the same manner as that predicted when the nonconservative 
forces are not present. It is assumed that in the damped system with random excitation, 
the natural frequency depends on the envelope of the response of the original nonlinear 
system. The present method thus depends on knowledge of the statistics of the response 
envelope of the original nonlinear system being studied. The probability density of the 
response envelope is known for several classes of nonlinear systems [10]. The random 
fluctuations m the natural frequency are then assumed to occur much more slowly than 
the fluctuations in the response. This leads to a simple formula for the response spectrum 
m the form of an integration over the probability density of the envelope of the response 
ot the original nonlinear oscillator. 

The present approach has been applied to a bilinear oscillator where the restoring 
force is assumed to be equal to the deflection ( stiffness = 1) for small deflections and is 
equal to twice the deflection (stiffness = 2) when the amplitude of the deflection is greater 
than unity. The damping is assumed to be linear viscous damping. This system is highly 
nonlinear when the random response spends a significant amount of time crossing over be- 
tween the two stiffness regimes. Numerical simulations of the random response have shown 
that this bilinear stiffness characteristic has a pronounced effect on the predicted power 
spectrum. It is hoped that if accurate spectral estimates can be obtained for this rather 
difficult system, then the approximate method may be applicable to a fairly broad class of 
oscillators. Comparisons of spectra obtained using the present approximate scheme with 
those obtained using numerical simulations are presented and show excellent agreement. 

As discussed above, the main purpose of the present study is to propose a description 
ol the influence of nonlinearities on the power spectra of random structures with random 
excitation. Although the basic assumptions of the approach are plausible, a detailed error 
analysis has yet to be conducted. The validity of the approach is investigated here through 
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a single example, the bilinear oscillator. A more detailed investigation of the influence of 
a number of approximations in the procedure will be performed in a future study. 


XI.2. Approximate Representation of the Nonlinear System 


Consider a nonlinear oscillator governed by 


x+u 0 2 (x + eg(x)) + cx = f(t ) (XI. 1) 

where f(t) is Gaussian white noise, g(x) is a nonlinear restoring force, u 0 is the natural 
frequency in the absence of nonlinear effects, e is a constant and c is a viscous damping 
coefficient. For the present investigation we consider only systems having conservative 
nonlinearities. The function g(x) does not depend on velocity, x. Nonlinear damping 
effects will be left for future studies. It is also assumed that g(x) is an odd function of x 
i.e. g(x) = -g(-x). 

To develop a representation of the nonlinear random response, consider the behavior of 
a conservative nonlinear system. It is well known that in the case of conservative nonlinear 
systems, where c and f(t) are zero, the solution of equation (XI.l), x(t), will consist of 
an oscillation having a period that depends on the amplitude of the motion. The exact 
solution for the period, T, corresponding to an oscillation amplitude, a, is given by 


T(a) = 4 



oj 0 2 (u + e<7(u))du] x ! 2 dx 


(XI.2) 


In the case of the damped system with random excitation described by equation (XI.l), 
if the damping coefficient, c, is sufficiently small the response may be considered to be a 
narrowband random process. The period of the oscillation cycles can be taken to be the 
time duration between occurrences of zero crossings with positive slope. The time required 
for one such cycle will depend on some measure of the oscillation amplitude. This follows 
from the fact that the stiffness of the nonlinear system, and hence, the effective natural 
frequency, depend on the response amplitude. 

One can also view the response in the phase plane. In the case of the conservative 
system where c and /(f) are zero, the amplitude of a given orbit in the phase plane deter- 
mines the oscillation frequency. If the damping and excitation are sufficiently small, during 
one cycle the orbit of the nonconservative system will consist of a small fluctuation about 
that of the conservative system. If the amplitudes of the conservative and nonconservative 
orbits are nearly the same, the time duration of the cycles will also be similar. 

During the forced, damped response of the system described by equation (XI.l), the 
amplitude of the motion will vary randomly and, based on the above discussion, we can 
also expect the oscillation period to vary accordingly. An approximate representation of 
the response of the system described by equation (XI.l) may then be obtained by idealizing 
it as a linear system having a natural frequency that depends on the response amplitude. 
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This amplitude dependent natural frequency is taken to be identical to that obtained for 
the conservative system, 

w(a) = 2tt /T(a), (XI.3) 

where T(a ) is given in equation (XI. 2). The oscillation frequency is thus regarded as an 
intrinsic property of the system; it does not depend explicitly on the excitation or the 
damping as long as the damping is sufficiently light. During the forced, damped response, 
the oscillation frequency depends only on the amplitude of the motion as in the case of 
the conservative system. The amplitude, a, is a random quantity which for the present 
study is taken to be equal to the envelope of the random response of the original system 
in equation (XI.l). If f(t) is Gaussian white noise, then the probability density of the 
envelope of the solution to equation (XI.l) may be determined [10]. Equation (XI.l) may 
then be approximated by 

x + u 2 (a)x + cx = f(t). (XI. 4) 


XI.3 Estimates of the Response Power Spectral Density 

Having a representation of the nonlinear system of equation (XI.l) in the form of a 
linear system with a randomly varying natural frequency as in equation (XI.4) allows the 
estimation of the response power spectral density. This may be obtained by taking the 
Fourier transform of the autocorrelation function of the stationary response, £'[r(t)x(t+r)l, 
where £[•] denotes the expected value. The solution of equation (XI.4) may be written in 
the form of a convolution integral, 


x (0 = J o h (t ~ T)f{r)dT, (XI.5) 

where h(t) is the impulse response which is assumed to depend on the randomly varying 
natural frequency, w(a), given in equation (XI.3). The estimation of the power spectral 
density proceeds m the usual manner (see for example reference [10]) with the exception 
that we must account for the random nature of the envelope, a. It is assumed that although 
a varies randomly, it varies much more slowly than x(<). An approximate expression for 
h{t) m equation (XI.5) may be obtained by solving 

h + u 2 (a)h + ch = 6(t), (XI.6) 

with initial conditions, A(0) = h( 0) = 0, and where £(•) is the Dirac delta function. This 
is equivalent to 

h + u 2 (a)h + ch = 0, (XI. 7) 

with initial conditions, h(0) = 0,A(0) = 1. If we assume that ^(a) varies much more 
slowly than h{t), the solution of equation (XI. 7) may be approximated by 


h(t) = 


e~ c< / 2 sm(V u> 2 (a)-(c /2) 2 ) 

>M«)-(c/2) 2 


(XI.8) 
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Given A(t), one can calculate the autocorrelation function as 


Rzx(t) - E[x(t)x(t+r)] = E[(J^ A(r , )/(<-r'Vr , )( j h(r")f{t + T-T")dT")]. (XI.9) 

Since the impulse response in equation (XI.8) is zero for negative values of t, the lower limits 
of integration may be set to -oo. Because we are interested in the stationary response, 
the value of t is set to oo. The upper limits of integration axe then set to oo. By assuming 
that the impulse response and the white noise excitation are uncorrelated, equation (XI.9) 
may be written as ' 


poo poo 

Rzz(t) = J ^ j ^ E[h(T')h(T n )]E[f(t - r')f(t + r - r")]dr , dr". (XI.10) 


It is assumed that f(t ) is stationary Gaussian white noise so that 

E[f(t - r')f(t + t- t")\ = 27 T$fS(T + t'~ t"), 


(XI.11) 


where S(-) is the Dirac delta function and is the power spectral density of the excitation. 
The integration over r" in equation (XI. 10) then gives 


Rxz( t ) — f E[h(T f )h(T + r^Jdr'. 

7 — 00 


(XI. 12) 


The power spectral density of x(t) is equal to the Fourier transform of R xx {r) , 
*.(“0 = ^J Rzz(T)e~ iwT dr 

r°° roo (XI. 13) 

= *// / E[h(r')h(T + T')]e- iu,T dT'dT, 

7 — oo 7 — oo 

where i is v/-L If the integrand is multiplied by e iu,T ' e~ iur> , equation (XI.13) becomes 

/ OO poo 

J + (XI.14) 

carrying out the integration first over t and then T r gives 

= $/£[5(w)ff*(«)], (XI. 15) 

where (•)* denotes the complex conjugate and, 

H{uj) = / h(r)e“ ,u,r dr. (XI.16) 

7 — oo 
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Substituting equation (XI.8) into equations (XI.15) and (XI. 16) gives the response 
spectrum in the form 

= ^(u> 2 (a) — a> 2 ^) 2 + (cu>) 2 ^’ (XL17 ) 

where in the evaluation of equation (XI. 16) we have again assumed that a; 2 (a) is slowly 
varying in comparison to x(i). 

The calculation of the expected value in equation (XI. 17) is possible if the probability 
density of the envelope response, a , is known. The exact expression for the probability 
density of the envelope of the solution to equation (XI. 1) is [10] 

P(a) = [u>l(x + e5(x))27r/(c„a;(a))]e“ V(a)/<r o, (XI.18) 

where c n is a normalization constant so that 

/ P(a)da = 1, (XI.19) 

J — OO 

and V(a) is the potential energy of the system corresponding to the response amplitude a, 

V(a) = u)%(a 2 /2 + e J g(x)dx). (XI.20) 

°o i s th e mean square response of equation (XI. 1) when e = 0, 

.2 _ */* 


cr 0 = 


aoi 


(XI.21) 


The approximate response spectrum is then given by 

^*0*0 = / T~ 2 (~\ ~ i \7 } — r^P(a)da. 

Jo (u> 2 (a) - w 2 ) 2 + (cw) 2 v 1 


(XI. 22) 


XI.4 Numerical Results for a Bilinear Oscillator 

Equation (XI.22) has been applied to estimate the response spectrum of a bilinear 
system having a nonlinear restoring force as shown in figure XI.l. The restoring force 
shown in figure XI.l is equal to w 2 (x + e ff (x)) in equation (XI.l). In the present calculations 
the system is assumed to have a linear force/deflection characteristic with unit slope when 
the amplitude of the response is less than unity. For greater response amplitudes, the 
stiffness of the system is assumed to be doubled as shown in the figure. The viscous 
damping coefficient, c, in equation (XI.l) is taken to be .005 . 

The results of applying equation (XI.22) for a range of input force spectrum levels, 
$/ are shown in figures XI.2 A) through D). The figures also show the estimated spectra 
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obtained by direct numerical simulation of the response of equation (XI. 1). This was 
computed by calculating the time domain response and estimating the spectrum by a 
ast Fourier Transform. The figures show that the power spectrum exhibits a broadened 
shape when the nonlinearity is significant. The effect of the nonlinear restoring force on 
the spectrum shape is depicted extremely well by the approximate method of equation 
(XI. 22). 

The primary discrepancy appears in figures XI.2 C) and XI.2 D) where the approxi- 
mate expression in equation (XI.22) predicts a significant peak a 1 rad/sec. Since neither 
solution is exact it is desirable to employ a third solution method to examine this effect. 
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XII Conclusions 


Several analytical methods have been developed to facilitate the prediction of the fatigue 
lives of structures that are subjected to acoustic loadings. In sections II through XI, time 
domain methods were presented for performing fatigue predictions based on knowledge of 
the power spectral density of the random strain response. It is suggested that procedures 
such as those presented here will lead to much more reliable fatigue predictions than 
methods that are commonly used. 

Because any fatigue prediction procedure relies on empirical data to characterize ma- 
terial fatigue properties, in section VII we present a technique for estimating the necessary 
constants using either narrowband random loads or realistic broadband loads as encoun- 
tered in service. One advantage of this procedure is that it enables one to identify fatigue 
constants from tests of complex structures with multiple resonant modes rather than sim- 
plified coupons. 

Another contribution of the present study has been to construct a method of account- 
ing for nonlinear response in random fatigue predictions. This can be helpful in situations 
where the in service loads are sufficiently intense to elicit nonlinear response or where non- 
linear effects are found to be significant in accelerated tests which employ artificially high 
excitation levels. The approximate procedure is applied to nonlinear beams and plates 
and excellent agreement is observed between results using the approximate method and 
using a detailed numerical simulation. Conventional numerical simulation methods are not 
practical for random fatigue predictions for nonlinear structures because the nonlinearities 
drastically complicate the calculations. It is found that the approximate method makes it 
practical to analyze the fatigue lives of highly complex nonlinear structures. This is be- 
cause the numerical effort required using the approximate procedure is roughly the same 
as for a linear structure. 

The final contribution of this effort has been to identify the effect of nonlinearities 
on the power spectral density of the random response. The power spectral density is 
commonly used in the analysis of random vibration. When nonlinear effects contribute to 
the response, the power spectrum is influenced in a manner that is far from obvious. An 
explanation of this effect is proposed and an approximate analytical method is presented to 
predict the power spectrum of nonlinear random systems. Excellent agreement is observed 
between the approximate method and numerical simulations. 
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